New class for computing corrections as a function of several variables (Chiara)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Feb 2009 01:27:30 +0000 (01:27 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Feb 2009 01:27:30 +0000 (01:27 +0000)
PWG3/CMake_libPWG3vertexingHF.txt
PWG3/PWG3vertexingHFLinkDef.h
PWG3/libPWG3vertexingHF.pkg
PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.C [new file with mode: 0644]
PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx [new file with mode: 0644]
PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.h [new file with mode: 0644]
PWG3/vertexingHF/ReadCFHeavyFlavourOutput.C [new file with mode: 0644]

index 7564d55..adfa995 100644 (file)
@@ -10,6 +10,7 @@ set(SRCS
        vertexingHF/AliAnalysisTaskSESelectHF.cxx
        vertexingHF/AliAnalysisTaskSECompareHF.cxx
        vertexingHF/AliCFHeavyFlavourTask.cxx
+       vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx
        vertexingHF/AliMultiDimVector.cxx
        vertexingHF/AliSignificanceCalculator.cxx
 )
index 3f85d2b..07bb437 100644 (file)
@@ -15,6 +15,7 @@
 #pragma link C++ class AliAnalysisTaskSESelectHF+;
 #pragma link C++ class AliAnalysisTaskSECompareHF+;
 #pragma link C++ class AliCFHeavyFlavourTask+;
+#pragma link C++ class AliCFHeavyFlavourTaskMultiVar+;
 #pragma link C++ class AliMultiDimVector+;
 #pragma link C++ class AliSignificanceCalculator+;
 
index 5101a10..2b7178f 100644 (file)
@@ -9,6 +9,7 @@ SRCS:=   vertexingHF/AliAODRecoDecayHF.cxx \
        vertexingHF/AliAnalysisTaskSESelectHF.cxx \
        vertexingHF/AliAnalysisTaskSECompareHF.cxx \
        vertexingHF/AliCFHeavyFlavourTask.cxx \
+       vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx \
        vertexingHF/AliMultiDimVector.cxx vertexingHF/AliSignificanceCalculator.cxx
 
 HDRS:= $(SRCS:.cxx=.h)
diff --git a/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.C b/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.C
new file mode 100644 (file)
index 0000000..8545a00
--- /dev/null
@@ -0,0 +1,294 @@
+//DEFINITION OF A FEW CONSTANTS
+const Double_t ymin  = -2.0 ;
+const Double_t ymax  =  2.0 ;
+const Double_t ptmin_0_4 =  0.0 ;
+const Double_t ptmax_0_4 =  4.0 ;
+const Double_t ptmin_4_8 =  4.0 ;
+const Double_t ptmax_4_8 =  8.0 ;
+const Double_t ptmin_8_10 =  8.0 ;
+const Double_t ptmax_8_10 =  10.0 ;
+const Double_t cosmin = -1.;
+const Double_t cosmax =  1.;
+const Double_t cTmin = 0;  // micron
+const Double_t cTmax = 500;  // micron
+const Int_t    mintrackrefsTPC = 2 ;
+const Int_t    mintrackrefsITS = 3 ;
+const Int_t    charge  = 1 ;
+const Int_t    PDG = 421; 
+const Int_t    minclustersTPC = 50 ;
+//----------------------------------------------------
+
+Bool_t AliCFHeavyFlavourTaskMultiVar()
+{
+       
+       TBenchmark benchmark;
+       benchmark.Start("AliCFHeavyFlavourTaskMultiVar");
+       
+       AliLog::SetGlobalDebugLevel(0);
+       
+       Load() ; //load the required libraries
+       
+
+       TChain * analysisChain ;
+       TChain * analysisChainFriend ;
+       
+       //here put your input data path
+       printf("\n\nRunning on local file, please check the path\n\n");
+       analysisChain = new TChain("aodTree");
+       analysisChainFriend = new TChain("aodTree");
+       for (Int_t i = 1; i<2; i++){
+               if (i != 4 && i != 5 && i != 6 && i != 10 && i != 11 && i != 14 && i != 15  && i != 16  && i != 18 && i != 20 && i != 21 && i != 24 && i != 26 && (i <30 || i > 43) && i!= 46 && i!=48 && i!=49 && i!=51 && i!=53 && i!=56 && i!=64 && i!=73 && i!=76 && i!=78 && i!=83 && i!=84 && i!=85 && i!=90 && i!=91 && i!=92){ 
+
+                       TString fileName = "/home/zampolli/HF/FromRenu/Events/180002/";
+                       TString aodHFFileName = "/home/zampolli/HF/FromRenu/Events/180002/";
+                       fileName+=Form("%i/AliAOD.root",i);
+                       printf(Form("file = %s \n", fileName.Data()));
+                       aodHFFileName+=Form("%i/AliAOD.VertexingHF.root",i);
+                       analysisChain->Add(fileName);
+                       analysisChainFriend->Add(aodHFFileName);
+               }
+       }
+                       
+       analysisChain->AddFriend(analysisChainFriend);
+       Info("AliCFHeavyFlavourTaskMultiVar",Form("CHAIN HAS %d ENTRIES",(Int_t)analysisChain->GetEntries()));
+       
+       //CONTAINER DEFINITION
+       Info("AliCFHeavyFlavourTaskMultiVar","SETUP CONTAINER");
+       //the sensitive variables (6 in this example), their indices
+       UInt_t ipt = 0;
+       UInt_t iy  = 1;
+       UInt_t icosThetaStar  = 2;
+       UInt_t ipTpi  = 3;
+       UInt_t ipTk  = 4;
+       UInt_t icT  = 5;
+       //Setting up the container grid... 
+       UInt_t nstep = 2 ; //number of selection steps MC 
+       const Int_t nvar   = 6 ; //number of variables on the grid:pt, y, cosThetaStar, pTpi, pTk, cT 
+       const Int_t nbin0_0_4  = 8 ; //bins in pt from 0 to 4 GeV
+       const Int_t nbin0_4_8  = 4 ; //bins in pt from 4 to 8 GeV
+       const Int_t nbin0_8_10  = 1 ; //bins in pt from 8 to 10 GeV
+       const Int_t nbin1  = 8 ; //bins in y
+       const Int_t nbin2  = 8 ; //bins in cosThetaStar 
+       const Int_t nbin3_0_4  = 8 ; //bins in ptPi from 0 to 4 GeV
+       const Int_t nbin3_4_8  = 4 ; //bins in ptPi from 4 to 8 GeV
+       const Int_t nbin3_8_10  = 1 ; //bins in ptPi from 8 to 10 GeV
+       const Int_t nbin4_0_4  = 8 ; //bins in ptKa from 0 to 4 GeV
+       const Int_t nbin4_4_8  = 4 ; //bins in ptKa from 4 to 8 GeV
+       const Int_t nbin4_8_10  = 1 ; //bins in ptKa from 8 to 10 GeV
+       const Int_t nbin5  = 24 ; //bins in cT
+       
+       //arrays for the number of bins in each dimension
+       Int_t iBin[nvar];
+       iBin[0]=nbin0_0_4+nbin0_4_8+nbin0_8_10;
+       iBin[1]=nbin1;
+       iBin[2]=nbin2;
+       iBin[3]=nbin3_0_4+nbin3_4_8+nbin3_8_10;
+       iBin[4]=nbin4_0_4+nbin4_4_8+nbin4_8_10;
+       iBin[5]=nbin5;
+       
+       //arrays for lower bounds :
+       Double_t *binLim0=new Double_t[iBin[0]+1];
+       Double_t *binLim1=new Double_t[iBin[1]+1];
+       Double_t *binLim2=new Double_t[iBin[2]+1];
+       Double_t *binLim3=new Double_t[iBin[3]+1];
+       Double_t *binLim4=new Double_t[iBin[4]+1];
+       Double_t *binLim5=new Double_t[iBin[5]+1];
+       
+       // checking limits
+       if (ptmax_0_4 != ptmin_4_8) {
+               Error("AliCFHeavyFlavourTaskMultiVar","max lim 1st range != min lim 2nd range, please check!");
+       }
+       if (ptmax_4_8 != ptmin_8_10) {
+               Error("AliCFHeavyFlavourTaskMultiVar","max lim 2nd range != min lim 3rd range, please check!");
+       }
+
+       // values for bin lower bounds
+       // pt
+       for(Int_t i=0; i<=nbin0_0_4; i++) binLim0[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin0_0_4*(Double_t)i ; 
+       if (binLim0[nbin0_0_4] != ptmin_4_8)  {
+               Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for pt - 1st range - differs from expected!\n");
+       }
+       for(Int_t i=0; i<=nbin0_4_8; i++) binLim0[i+nbin0_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin0_4_8*(Double_t)i ; 
+       if (binLim0[nbin0_0_4+nbin0_4_8] != ptmin_8_10)  {
+               Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for pt - 2nd range - differs from expected!\n");
+       }
+       for(Int_t i=0; i<=nbin0_8_10; i++) binLim0[i+nbin0_0_4+nbin0_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin0_8_10*(Double_t)i ; 
+
+       // y
+       for(Int_t i=0; i<=nbin1; i++) binLim1[i]=(Double_t)ymin  + (ymax-ymin)  /nbin1*(Double_t)i ;
+
+       // cosThetaStar
+       for(Int_t i=0; i<=nbin2; i++) binLim2[i]=(Double_t)cosmin  + (cosmax-cosmin)  /nbin2*(Double_t)i ;
+
+       // ptPi
+       for(Int_t i=0; i<=nbin3_0_4; i++) binLim3[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin3_0_4*(Double_t)i ; 
+       if (binLim3[nbin3_0_4] != ptmin_4_8)  {
+               Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for ptPi - 1st range - differs from expected!");
+       }
+       for(Int_t i=0; i<=nbin3_4_8; i++) binLim3[i+nbin3_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin3_4_8*(Double_t)i ; 
+       if (binLim3[nbin3_0_4+nbin3_4_8] != ptmin_8_10)  {
+               Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for ptPi - 2nd range - differs from expected!\n");
+       }
+       for(Int_t i=0; i<=nbin3_8_10; i++) binLim3[i+nbin3_0_4+nbin3_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin3_8_10*(Double_t)i ; 
+
+       // ptKa
+       for(Int_t i=0; i<=nbin4_0_4; i++) binLim4[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin4_0_4*(Double_t)i ; 
+       if (binLim4[nbin4_0_4] != ptmin_4_8)  {
+               Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for ptKa - 1st range - differs from expected!");
+       }
+       for(Int_t i=0; i<=nbin4_4_8; i++) binLim4[i+nbin4_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin4_4_8*(Double_t)i ; 
+       if (binLim4[nbin4_0_4+nbin4_4_8] != ptmin_8_10)  {
+               Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for ptKa - 2nd range - differs from expected!\n");
+       }
+       for(Int_t i=0; i<=nbin4_8_10; i++) binLim4[i+nbin4_0_4+nbin4_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin4_8_10*(Double_t)i ; 
+
+       // cT
+       for(Int_t i=0; i<=nbin5; i++) binLim5[i]=(Double_t)cTmin  + (cTmax-cTmin)  /nbin5*(Double_t)i ;
+
+       Info("AliCFHeavyFlavourTaskMultiVar","Printing lower limits for bins in pt");
+       for (Int_t i =0; i<= iBin[0]; i++){
+               Info("AliCFHeavyFlavourTaskMultiVar",Form("i-th bin, lower limit = %f", binLim0[i]));
+       }
+       Info("Printing lower limits for bins in ptPi");
+       for (Int_t i =0; i<= iBin[3]; i++){
+               Info("AliCFHeavyFlavourTaskMultiVar",Form("i-th bin, lower limit = %f", binLim3[i]));
+       }
+       Info("Printing lower limits for bins in ptKa");
+       for (Int_t i =0; i<= iBin[4]; i++){
+               Info("AliCFHeavyFlavourTaskMultiVar",Form("i-th bin, lower limit = %f", binLim4[i]));
+       }
+
+       //one "container" for MC
+       AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin);
+       //setting the bin limits
+       container -> SetBinLimits(ipt,binLim0);
+       container -> SetBinLimits(iy,binLim1);
+       container -> SetBinLimits(icosThetaStar,binLim2);
+       container -> SetBinLimits(ipTpi,binLim3);
+       container -> SetBinLimits(ipTk,binLim4);
+       container -> SetBinLimits(icT,binLim5);
+       
+       //CREATE THE  CUTS -----------------------------------------------
+       
+       // Gen-Level kinematic cuts
+       AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
+       //   mcKineCuts->SetPtRange(ptmin,ptmax);
+       //   mcKineCuts->SetRapidityRange(ymin,ymax);
+       //   mcKineCuts->SetChargeMC(charge);
+       
+       //Particle-Level cuts:  
+       AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
+       //mcGenCuts->SetRequireIsPrimary();
+       mcGenCuts->SetRequirePdgCode(PDG, kTRUE);  // kTRUE set in order to include D0_bar
+       mcGenCuts->SetAODMC(1); //special flag for reading MC in AOD tree (important)
+       
+       // Rec-Level kinematic cuts
+       AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");
+       //   recKineCuts->SetPtRange(ptmin,ptmax);
+       //   recKineCuts->SetRapidityRange(ymin,ymax);
+       //   recKineCuts->SetChargeRec(charge);
+       
+       AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
+       //recQualityCuts->SetStatus(AliESDtrack::kITSrefit);
+       
+       AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
+       //recIsPrimaryCuts->SetAODType(AliAODTrack::kPrimary);
+       
+       printf("CREATE MC KINE CUTS\n");
+       TObjArray* mcList = new TObjArray(0) ;
+       mcList->AddLast(mcKineCuts);
+       mcList->AddLast(mcGenCuts);
+       
+       printf("CREATE RECONSTRUCTION CUTS\n");
+       TObjArray* recList = new TObjArray(0) ;
+       recList->AddLast(recKineCuts);
+       recList->AddLast(recQualityCuts);
+       recList->AddLast(recIsPrimaryCuts);
+       
+       //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
+       printf("CREATE INTERFACE AND CUTS\n");
+       AliCFManager* man = new AliCFManager() ;
+       man->SetParticleContainer     (container);
+       man->SetParticleCutsList(0 , mcList); // MC
+       man->SetParticleCutsList(1 , recList); // AOD
+       
+       //CREATE THE TASK
+       printf("CREATE TASK\n");
+       // create the task
+       AliCFHeavyFlavourTaskMultiVar *task = new AliCFHeavyFlavourTaskMultiVar("AliCFHeavyFlavourTaskMultiVar");
+       task->SetFillFromGenerated(kFALSE);
+       task->SetCFManager(man); //here is set the CF manager
+       
+       //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS
+       printf("CREATE ANALYSIS MANAGER\n");
+       // Make the analysis manager
+       AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+       
+       mgr->SetAnalysisType(AliAnalysisManager::kLocalAnalysis);
+       
+       AliAODInputHandler* dataHandler = new AliAODInputHandler();
+       mgr->SetInputEventHandler(dataHandler);
+       
+       // Create and connect containers for input/output
+       
+       // ------ input data ------
+       AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();
+       
+       // ----- output data -----
+       
+       //slot 0 : default output tree (by default handled by AliAnalysisTaskSE)
+       AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("ctree0", TTree::Class(),AliAnalysisManager::kOutputContainer,"output.root");
+       
+       //now comes user's output objects :
+       
+       // output TH1I for event counting
+       AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,"output.root");
+       // output Correction Framework Container (for acceptance & efficiency calculations)
+       AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,"output.root");
+       
+       //cinput0->SetData(analysisChain);
+       
+       mgr->AddTask(task);
+       
+       mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+       mgr->ConnectOutput(task,0,coutput0);
+       mgr->ConnectOutput(task,1,coutput1);
+       mgr->ConnectOutput(task,2,coutput2);
+       
+       /*
+       mgr->ConnectOutput(task,0,mgr->GetCommonOutputContainer());
+       mgr->ConnectOutput(task,1,mgr->GetCommonOutputContainer());
+       mgr->ConnectOutput(task,2,mgr->GetCommonOutputContainer());
+       */
+
+       printf("READY TO RUN\n");
+       //RUN !!!
+       if (mgr->InitAnalysis()) {
+               mgr->PrintStatus();
+               mgr->StartAnalysis("local",analysisChain);
+       }
+       
+       benchmark.Stop("AliCFHeavyFlavourTaskMultiVar");
+       benchmark.Show("AliCFHeavyFlavourTaskMultiVar");
+       
+       return kTRUE ;
+}
+
+void Load() {
+       
+       //load the required aliroot libraries
+       gSystem->Load("libANALYSIS") ;
+       gSystem->Load("libANALYSISalice") ;
+       gSystem->Load("libCORRFW.so") ;
+       gSystem->Load("libPWG3base.so");
+       gSystem->Load("libPWG3vertexingHF.so");
+
+       gSystem->Load("libTree.so");
+       gSystem->Load("libGeom.so");
+       gSystem->Load("libPhysics.so");
+       gSystem->Load("libVMC.so");
+       gSystem->Load("libSTEERBase.so");
+       gSystem->Load("libESD.so");
+       gSystem->Load("libAOD.so"); 
+       
+}
diff --git a/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx b/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx
new file mode 100644 (file)
index 0000000..a1d01d8
--- /dev/null
@@ -0,0 +1,806 @@
+/**************************************************************************
+ * 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 for HF corrections as a function of many variables
+//
+// Example of task running on AliEn
+// which provides standard way of calculating acceptance and efficiency
+// between different steps of the procedure.
+// The ouptut of the task is a AliCFContainer from which the efficiencies
+// can be calculated
+//-----------------------------------------------------------------------
+// Author : C. Zampolli, CERN, on the basis of R. Vernet's example
+//-----------------------------------------------------------------------
+
+#include <TCanvas.h>
+#include <TParticle.h>
+#include <TH1I.h>
+#include <TStyle.h>
+
+#include "AliCFHeavyFlavourTaskMultiVar.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliCFManager.h"
+#include "AliCFContainer.h"
+#include "AliLog.h"
+#include "AliAODEvent.h"
+#include "AliAODRecoDecay.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODMCParticle.h"
+
+//__________________________________________________________________________
+AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar() :
+       AliAnalysisTaskSE(),
+       fPDG(0),
+       fCFManager(0x0),
+       fHistEventsProcessed(0x0),
+       fCountMC(0),
+       fEvents(0),
+       fFillFromGenerated(kFALSE),
+       fSkipped(0)
+{
+       //
+       //Default ctor
+       //
+}
+//___________________________________________________________________________
+AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar(const Char_t* name) :
+       AliAnalysisTaskSE(name),
+       fPDG(0),
+       fCFManager(0x0),
+       fHistEventsProcessed(0x0),
+       fCountMC(0),
+       fEvents(0),
+       fFillFromGenerated(kFALSE),
+       fSkipped(0)
+{
+       //
+       // Constructor. Initialization of Inputs and Outputs
+       //
+       Info("AliCFHeavyFlavourTaskMultiVar","Calling Constructor");
+       /*
+         DefineInput(0) and DefineOutput(0)
+         are taken care of by AliAnalysisTaskSE constructor
+       */
+       DefineOutput(1,TH1I::Class());
+       DefineOutput(2,AliCFContainer::Class());
+}
+
+//___________________________________________________________________________
+AliCFHeavyFlavourTaskMultiVar& AliCFHeavyFlavourTaskMultiVar::operator=(const AliCFHeavyFlavourTaskMultiVar& c) 
+{
+       //
+       // Assignment operator
+       //
+       if (this!=&c) {
+               AliAnalysisTaskSE::operator=(c) ;
+               fPDG      = c.fPDG;
+               fCFManager  = c.fCFManager;
+               fHistEventsProcessed = c.fHistEventsProcessed;
+       }
+       return *this;
+}
+
+//___________________________________________________________________________
+AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar(const AliCFHeavyFlavourTaskMultiVar& c) :
+       AliAnalysisTaskSE(c),
+       fPDG(c.fPDG),
+       fCFManager(c.fCFManager),
+       fHistEventsProcessed(c.fHistEventsProcessed),
+       fCountMC(c.fCountMC),
+       fEvents(c.fEvents),
+       fFillFromGenerated(c.fFillFromGenerated),
+       fSkipped(c.fSkipped)
+{
+       //
+       // Copy Constructor
+       //
+}
+
+//___________________________________________________________________________
+AliCFHeavyFlavourTaskMultiVar::~AliCFHeavyFlavourTaskMultiVar() {
+       //
+       //destructor
+       //
+       Info("~AliCFHeavyFlavourTaskMultiVar","Calling Destructor");
+       if (fCFManager)           delete fCFManager ;
+       if (fHistEventsProcessed) delete fHistEventsProcessed ;
+}
+
+//_________________________________________________
+void AliCFHeavyFlavourTaskMultiVar::UserExec(Option_t *)
+{
+       //
+       // Main loop function
+       //
+       
+       if (!fInputEvent) {
+               Error("UserExec","NO EVENT FOUND!");
+               return;
+       }
+       
+       fEvents++;
+       if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
+       AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
+       fCFManager->SetEventInfo(aodEvent);
+       
+       // MC-event selection
+       Double_t containerInput[6] ;
+        
+       //loop on the MC event
+       
+       TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+       if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
+       Int_t icountMC = 0;
+       
+       for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { 
+               AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
+               if (!mcPart) {
+                       AliWarning("Particle not found in tree, skipping"); 
+                       continue;
+               } 
+               
+               // check the MC-level cuts
+               if (!fCFManager->CheckParticleCuts(0, mcPart)) continue;  // 0 stands for MC level
+               
+               // fill the container for Gen-level selection
+               Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.};
+               if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
+                       containerInput[0] = vectorMC[0];
+                       containerInput[1] = vectorMC[1] ;
+                       containerInput[2] = vectorMC[2] ;
+                       containerInput[3] = vectorMC[3] ;
+                       containerInput[4] = vectorMC[4] ;
+                       containerInput[5] = vectorMC[5] ;  // in micron
+                       fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
+                       icountMC++;
+               }
+               else {
+                       AliDebug(3,"Problems in filling the container");
+                       continue;
+               }
+       }    
+       
+       AliDebug(2, Form("Found %i MC particles that are D0!!",icountMC));
+       AliDebug(2, Form("aodEvent = %p ",aodEvent));
+
+       // Now go to rec level
+       fCountMC += icountMC;
+       // load heavy flavour vertices
+
+       TClonesArray *arrayVerticesHF = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi"));  
+       if (!arrayVerticesHF) AliError("Could not find array of HF vertices");
+       AliDebug(2, Form("Found %d vertices",arrayVerticesHF->GetEntriesFast()));
+       
+       for (Int_t iVertex = 0; iVertex<arrayVerticesHF->GetEntriesFast(); iVertex++) {
+               
+               AliAODRecoDecayHF2Prong* vtx = (AliAODRecoDecayHF2Prong*)arrayVerticesHF->At(iVertex);
+               
+               // cuts can't be applied to RecoDecays particles
+               // if (!fCFManager->CheckParticleCuts(1  , vtx)) continue;  // 1 stands for AOD level
+               
+               // find associated MC particle
+               Int_t mcLabel = IsMcVtx(vtx) ;
+               if (mcLabel == -1) 
+                       {
+                               AliDebug(2,"No MC particle found");
+                       }
+               else {
+                       AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
+                       if (!mcVtxHF) {
+                               AliWarning("Could not find associated MC in AOD MC tree");
+                               continue;
+                       }
+                                       
+                       // check if associated MC v0 passes the cuts
+                       if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) {  // 0 stands for MC
+                               AliDebug(2, "Skipping the particles due to cuts");
+                               continue; 
+                       }
+
+                       // fill the container
+                       // either with reconstructed values....
+                       Double_t pt = vtx->Pt();
+                       Double_t rapidity = vtx->YD0();
+                       
+                       Double_t cosThetaStar = 9999.;
+                       Double_t pTpi = 0.;
+                       Double_t pTK = 0.;
+                       Int_t pdgCode = mcVtxHF->GetPdgCode();
+                       if (pdgCode > 0){
+                               cosThetaStar = vtx->CosThetaStarD0();
+                               pTpi = vtx->PtProng(0);
+                               pTK = vtx->PtProng(1);
+                       }
+                       else {
+                               cosThetaStar = vtx->CosThetaStarD0bar();
+                               pTpi = vtx->PtProng(1);
+                               pTK = vtx->PtProng(0);
+                       }
+
+                       Double_t cT = vtx->CtD0();
+                       AliDebug(2, Form("cT from reconstructed vertex = %f (micron)",cT*1E4));
+                       AliDebug(2, Form("pdg code from MC = %d",TMath::Abs(mcVtxHF->GetPdgCode())));
+
+                       // ... or with generated values
+
+                       if (!fFillFromGenerated){
+                               containerInput[0] = pt;
+                               containerInput[1] = rapidity;
+                               containerInput[2] = cosThetaStar;
+                               containerInput[3] = pTpi;
+                               containerInput[4] = pTK;
+                               containerInput[5] = cT*1.E4;  // in micron
+                       }
+                       else {
+                               Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.};
+                               if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
+                                       containerInput[0] = vectorMC[0];
+                                       containerInput[1] = vectorMC[1] ;
+                                       containerInput[2] = vectorMC[2] ;
+                                       containerInput[3] = vectorMC[3] ;
+                                       containerInput[4] = vectorMC[4] ;
+                                       containerInput[5] = vectorMC[5] ;  // in micron
+                               }
+                               else {
+                                       AliDebug(3,"Problems in filling the container");
+                                       continue;
+                               }
+                       }
+                       AliDebug(2, Form("Filling the container with pt = %f, rapidity = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, cT = %f", containerInput[0], containerInput[1], containerInput[2], containerInput[3], containerInput[4], containerInput[5]));
+                       fCFManager->GetParticleContainer()->Fill(containerInput,1) ;   
+               }
+               /*
+               // for debugging only, filling container with dummy values
+               Double_t pt = 1.5;
+               Double_t rapidity = 0.5;
+               Double_t cosThetaStar = 0.7;
+               Double_t pTpi = 2.5;
+               Double_t pTK = 1;
+               Double_t cT = 3;
+               containerInput[0] = pt;
+               containerInput[1] = rapidity;
+               containerInput[2] = cosThetaStar;
+               containerInput[3] = pTpi;
+               containerInput[4] = pTK;
+               containerInput[5] = cT;
+               fCFManager->GetParticleContainer()->Fill(containerInput,1) ;   
+               */
+
+       }
+       
+       fHistEventsProcessed->Fill(0);
+       /* PostData(0) is taken care of by AliAnalysisTaskSE */
+       PostData(1,fHistEventsProcessed) ;
+       PostData(2,fCFManager->GetParticleContainer()) ;
+}
+
+
+//___________________________________________________________________________
+void AliCFHeavyFlavourTaskMultiVar::Terminate(Option_t*)
+{
+       // The Terminate() function is the last function to be called during
+       // a query. It always runs on the client, it can be used to present
+       // the results graphically or save the results to file.
+       
+       Info("Terminate","");
+       AliAnalysisTaskSE::Terminate();
+       
+       AliInfo(Form("Found %i MC particles that are D0 in MC in %d events",fCountMC,fEvents));
+       AliInfo(Form("Found %i reco D0 that are not decaying in K+pi in %d events",fSkipped,fEvents));
+       
+       // draw some example plots....
+       
+       AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
+       
+       // projecting the containers to obtain histograms
+       // first argument = variable, second argument = step
+       TH1D* h00 =   cont->ShowProjection(0,0) ;   // pt
+       TH1D* h10 =   cont->ShowProjection(1,0) ;   // rapidity
+       TH1D* h20 =   cont->ShowProjection(2,0) ;   // cosThetaStar
+       TH1D* h30 =   cont->ShowProjection(3,0) ;   // pTpi
+       TH1D* h40 =   cont->ShowProjection(4,0) ;   // pTK
+       TH1D* h50 =   cont->ShowProjection(5,0) ;   // cT
+       
+       TH1D* h01 =   cont->ShowProjection(0,1) ;   // pt
+       TH1D* h11 =   cont->ShowProjection(1,1) ;   // rapidity
+       TH1D* h21 =   cont->ShowProjection(2,1) ;   // cosThetaStar
+       TH1D* h31 =   cont->ShowProjection(3,1) ;   // pTpi
+       TH1D* h41 =   cont->ShowProjection(4,1) ;   // pTK
+       TH1D* h51 =   cont->ShowProjection(5,1) ;   // cT
+       
+       h00->SetTitle("pT_D0 (GeV/c)");
+       h10->SetTitle("rapidity");
+       h20->SetTitle("cosThetaStar");
+       h30->SetTitle("pT_pi (GeV/c)");
+       h40->SetTitle("pT_K (Gev/c)");
+       h50->SetTitle("cT (#mum)");
+
+       h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
+       h10->GetXaxis()->SetTitle("rapidity");
+       h20->GetXaxis()->SetTitle("cosThetaStar");
+       h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
+       h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
+       h50->GetXaxis()->SetTitle("cT (#mum)");
+
+       h01->SetTitle("pT_D0 (GeV/c)");
+       h11->SetTitle("rapidity");
+       h21->SetTitle("cosThetaStar");
+       h31->SetTitle("pT_pi (GeV/c)");
+       h41->SetTitle("pT_K (Gev/c)");
+       h51->SetTitle("cT (#mum)");
+
+       h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
+       h11->GetXaxis()->SetTitle("rapidity");
+       h21->GetXaxis()->SetTitle("cosThetaStar");
+       h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
+       h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
+       h51->GetXaxis()->SetTitle("cT (#mum)");
+
+       Double_t max0 = h00->GetMaximum();
+       Double_t max1 = h10->GetMaximum();
+       Double_t max2 = h20->GetMaximum();
+       Double_t max3 = h30->GetMaximum();
+       Double_t max4 = h40->GetMaximum();
+       Double_t max5 = h50->GetMaximum();
+       
+       h00->GetYaxis()->SetRangeUser(0,max0*1.2);
+       h10->GetYaxis()->SetRangeUser(0,max1*1.2);
+       h20->GetYaxis()->SetRangeUser(0,max2*1.2);
+       h30->GetYaxis()->SetRangeUser(0,max3*1.2);
+       h40->GetYaxis()->SetRangeUser(0,max4*1.2);
+       h50->GetYaxis()->SetRangeUser(0,max5*1.2);
+       
+       h01->GetYaxis()->SetRangeUser(0,max0*1.2);
+       h11->GetYaxis()->SetRangeUser(0,max1*1.2);
+       h21->GetYaxis()->SetRangeUser(0,max2*1.2);
+       h31->GetYaxis()->SetRangeUser(0,max3*1.2);
+       h41->GetYaxis()->SetRangeUser(0,max4*1.2);
+       h51->GetYaxis()->SetRangeUser(0,max5*1.2);
+       
+       h00->SetMarkerStyle(20);
+       h10->SetMarkerStyle(24);
+       h20->SetMarkerStyle(21);
+       h30->SetMarkerStyle(25);
+       h40->SetMarkerStyle(27);
+       h50->SetMarkerStyle(28);
+
+       h00->SetMarkerColor(2);
+       h10->SetMarkerColor(2);
+       h20->SetMarkerColor(2);
+       h30->SetMarkerColor(2);
+       h40->SetMarkerColor(2);
+       h50->SetMarkerColor(2);
+
+       h01->SetMarkerStyle(20) ;
+       h11->SetMarkerStyle(24) ;
+       h21->SetMarkerStyle(21) ;
+       h31->SetMarkerStyle(25) ;
+       h41->SetMarkerStyle(27) ;
+       h51->SetMarkerStyle(28) ;
+
+       h01->SetMarkerColor(4);
+       h11->SetMarkerColor(4);
+       h21->SetMarkerColor(4);
+       h31->SetMarkerColor(4);
+       h41->SetMarkerColor(4);
+       h51->SetMarkerColor(4);
+       
+       gStyle->SetCanvasColor(0);
+       gStyle->SetFrameFillColor(0);
+       gStyle->SetTitleFillColor(0);
+       gStyle->SetStatColor(0);
+
+       // drawing in 2 separate canvas for a matter of clearity
+       TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",800,1600);
+       TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",800,1600);
+       c1->Divide(2,3);
+       c2->Divide(2,3);
+       
+       c1->cd(1);
+       h00->Draw("p");
+       c1->cd(2);
+       h01->Draw("p");
+       c1->cd(3);
+       h10->Draw("p");
+       c1->cd(4);
+       h11->Draw("p");
+       c1->cd(5);
+       h20->Draw("p");
+       c1->cd(6);
+       h21->Draw("p");
+       c1->cd();
+
+       c2->cd(1);
+       h30->Draw("p");
+       c2->cd(2);
+       h31->Draw("p");
+       c2->cd(3);
+       h40->Draw("p");
+       c2->cd(4);
+       h41->Draw("p");
+       c2->cd(5);
+       h50->Draw("p");
+       c2->cd(6);
+       h51->Draw("p");
+       c2->cd();
+
+       c1->SaveAs("Variables/pT_rapidity_cosThetaStar.eps");
+       c2->SaveAs("Variables/pTpi_pTK_cT.eps");
+       
+}
+
+//___________________________________________________________________________
+void AliCFHeavyFlavourTaskMultiVar::UserCreateOutputObjects() {
+       //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
+       //TO BE SET BEFORE THE EXECUTION OF THE TASK
+       //
+       Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
+       
+       //slot #1
+       OpenFile(1);
+       fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
+}
+
+//___________________________________________________________________________
+Int_t AliCFHeavyFlavourTaskMultiVar::GetVtxLabel(Int_t* labels) {
+       //
+       // returns the label of the vertex, given the labels of the 2 daughter tracks
+       // returns -1 if the V0 is fake
+       //
+       
+       Int_t labD0daugh0=-1; 
+       Int_t labD0daugh1=-1;
+       Int_t labD0=-1;
+       Int_t labMother = -1;
+       Int_t pdgMother = -1;
+       
+       AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
+       TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+       if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
+
+       AliAODMCParticle *part0 = (AliAODMCParticle*)mcArray->At(labels[0]);
+       if(!part0) { 
+               printf("no MC particle\n");
+               return -1;
+       }
+
+       while(part0->GetMother()>=0) {
+               labMother=part0->GetMother();
+               part0 = (AliAODMCParticle*)mcArray->At(labMother);
+               if(!part0) {
+                       printf("no MC mother particle\n");
+               }
+               pdgMother = TMath::Abs(part0->GetPdgCode());
+               AliDebug(2,Form("pdgMother from particle 0 = %d",pdgMother));
+               if(pdgMother==421) {
+                       AliDebug(2,Form("found one particle (first daughter) coming from a D0"));
+                       labD0daugh0=labMother;
+               }
+       }
+       
+       AliAODMCParticle *part1 = (AliAODMCParticle*)mcArray->At(labels[1]);
+       if(!part1) {
+               printf("no MC particle\n");
+               return -1;
+       }
+
+       while(part1->GetMother()>=0) {
+               labMother=part1->GetMother();
+               part1 = (AliAODMCParticle*)mcArray->At(labMother);
+               if(!part1) {
+                       printf("no MC mother particle\n");
+               }
+               pdgMother = TMath::Abs(part1->GetPdgCode());
+               AliDebug(2,Form("pdgMother from particle 1 = %d",pdgMother));
+               if(pdgMother==421) {
+                       AliDebug(2,Form("found one particle (second daughter) coming from a D0"));
+                       labD0daugh1=labMother;
+               }
+       }
+       if ((labD0daugh0 >= 0) && (labD0daugh1 >= 0) && (labD0daugh0 == labD0daugh1)){
+               //AliDebug(2, Form("found particles from the same D0 mother"));
+               labD0 = labD0daugh0;
+               // check that the D0 decays in 2 prongs
+               AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(labD0);
+               if (TMath::Abs(part->GetDaughter(1)-part->GetDaughter(0))==1) {
+                      return labD0;
+               }
+               else {
+                       fSkipped+=1;
+               }
+       }
+       
+       return -1;
+
+}
+
+//___________________________________________________________________________
+Int_t AliCFHeavyFlavourTaskMultiVar::IsMcVtx(AliAODRecoDecayHF2Prong* vtx) {
+       //
+       // check if the passed vertex is associated to a MC one, 
+       //   and returns the corresponding geant label.
+       // returns -1 if it is fake (i.e. label<0).
+       //
+
+       
+       AliAODTrack *trk0 = (AliAODTrack*)vtx->GetDaughter(0);
+       AliAODTrack *trk1 = (AliAODTrack*)vtx->GetDaughter(1);
+
+       if (!trk0 || ! trk1) {
+               AliDebug(2, "problems with the tracks while looking for mother");
+               //if (!trk0) AliDebug(2, "track 0 gives problems");
+               //if (!trk1) AliDebug(2, "track 1 gives problems");
+               return -1;
+       }
+
+       Int_t labels[2] = {-1,-1};
+       labels[0] = trk0->GetLabel();
+       labels[1] = trk1->GetLabel();
+
+       if (labels[0] < 0 || labels[1] < 0) {
+               AliDebug(2, "problems with the labels of the tracks while looking for mother");
+               return -1;
+       }
+       
+       return GetVtxLabel(&labels[0]) ;
+}
+//___________________________________________________________________________
+Double_t AliCFHeavyFlavourTaskMultiVar::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
+
+       //
+       // to calculate cos(ThetaStar) for generated particle
+       // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively 
+       // (see where the function is called)
+       //
+
+       Int_t pdgvtx = mcPart->GetPdgCode();
+       /*      if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
+               Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
+               Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
+               AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
+               AliDebug(2,"This is a D0");
+               AliAODMCParticle* mcPartdummy = mcPartDaughter0;
+               mcPartDaughter0 = mcPartDaughter1;
+               mcPartDaughter1 = mcPartdummy;
+       } 
+       else{
+               AliInfo("D0bar");
+       }
+       */
+       Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
+       Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
+       if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
+               AliDebug(2,"D0");
+       }
+       else{
+               AliDebug(2,"D0bar");
+       }
+       if (pdgprong0 == 211){
+               AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
+               AliAODMCParticle* mcPartdummy = mcPartDaughter0;
+               mcPartDaughter0 = mcPartDaughter1;
+               mcPartDaughter1 = mcPartdummy;
+               pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
+               pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
+       } 
+
+       AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
+       Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
+       Double_t massp[2];
+       massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
+       massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
+
+       Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
+
+       Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
+       Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
+       Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
+       Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
+       Double_t e =  TMath::Sqrt(massvtx*massvtx+p*p);
+       Double_t beta = p/e;
+       Double_t gamma = e/massvtx;
+       TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
+       TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
+
+       Double_t qlprong = mom.Dot(momTot)/momTot.Mag();  // analog to AliAODRecoDecay::QlProng(0)
+       
+       AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
+       Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
+       AliDebug(2,Form("cts = %f", cts));
+       return cts;
+}
+//___________________________________________________________________________
+Double_t AliCFHeavyFlavourTaskMultiVar::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
+
+       //
+       // to calculate cT for generated particle
+       //
+
+       Double_t xmcPart[3] = {0,0,0};
+       Double_t xdaughter0[3] = {0,0,0};
+       Double_t xdaughter1[3] = {0,0,0};
+       mcPart->XvYvZv(xmcPart);  // cm
+       mcPartDaughter0->XvYvZv(xdaughter0);  // cm
+       mcPartDaughter1->XvYvZv(xdaughter1);  //cm
+       Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
+                                        xmcPart[1]*xmcPart[1]+
+                                        xmcPart[2]*xmcPart[2]);
+       Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
+                                               xdaughter0[1]*xdaughter0[1]+
+                                               xdaughter0[2]*xdaughter0[2]);
+       Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
+                                               xdaughter1[1]*xdaughter1[1]+
+                                               xdaughter1[2]*xdaughter1[2]);
+       AliDebug(2, Form("D0:        x = %f, y = %f, z = %f, production vertex distance = %f (cm), %f (micron)", xmcPart[0], xmcPart[1], xmcPart[2], prodVtxD0, prodVtxD0*1.E4));
+       AliDebug(2, Form("Daughter0: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter0[0], xdaughter0[1], xdaughter0[2], prodVtxDaughter0, prodVtxDaughter0*1E4));
+       AliDebug(2, Form("Daughter1: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter1[0], xdaughter1[1], xdaughter1[2], prodVtxDaughter1, prodVtxDaughter1*1.E4));
+
+       Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
+                                    (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
+                                    (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
+
+       Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
+                                    (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
+                                    (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
+
+       if (cT0 != cT1) {
+               AliWarning("cT from daughter 0 different from cT from daughter 1! Using cT from daughter 0, but PLEASE, CHECK....");
+       }
+
+       // calculating cT from cT0
+
+       Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
+       Double_t p = mcPart-> P();
+       Double_t cT = cT0*mass/p;
+       AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4)); 
+       AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4)); 
+       AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
+       return cT;
+}
+//________________________________________________________________________________
+Bool_t AliCFHeavyFlavourTaskMultiVar::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
+
+       // 
+       // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
+       //
+
+       Bool_t isOk = kFALSE;
+
+       // check whether the D0 decays in pi+K
+       // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       // could use a cut in the CF, but not clear how to define a TDecayChannel
+       // implemented for the time being as a cut on the number of daughters - see later when 
+       // getting the daughters
+
+       // getting the daughters
+       Int_t daughter0 = mcPart->GetDaughter(0);
+       Int_t daughter1 = mcPart->GetDaughter(1);
+       AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
+       if (daughter0 == 0 || daughter1 == 0) {
+               AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
+               return isOk;  
+       }
+       if (TMath::Abs(daughter1 - daughter0 != 1)) {
+               AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
+               return isOk;  
+       }
+       AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
+       AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
+       if (!mcPartDaughter0 || !mcPartDaughter1) {
+               AliWarning("At least one Daughter Particle not found in tree, skipping"); 
+               return isOk;  
+       }
+       
+       Double_t vtx1[3] = {0,0,0};   // primary vertex         
+       Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
+       Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
+       mcPart->XvYvZv(vtx1);  // cm
+       // getting vertex from daughters
+       mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
+       mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
+       if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
+               AliError("Daughters have different secondary vertex, skipping the track");
+               return isOk;
+       }
+       Int_t nprongs = 2;
+       Short_t charge = 0;
+       // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
+       AliAODMCParticle* positiveDaugh = mcPartDaughter0;
+       AliAODMCParticle* negativeDaugh = mcPartDaughter1;
+       if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
+               // inverting in case the positive daughter is the second one
+               positiveDaugh = mcPartDaughter1;
+               negativeDaugh = mcPartDaughter0;
+       }
+       // getting the momentum from the daughters
+       Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};            
+       Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};            
+       Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
+
+       Double_t d0[2] = {0.,0.};               
+
+       AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
+
+       Double_t cosThetaStar = 0.;
+       Double_t cosThetaStarD0 = 0.;
+       Double_t cosThetaStarD0bar = 0.;
+       cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
+       cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
+       if (mcPart->GetPdgCode() == 421){  // D0
+               AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
+               cosThetaStar = cosThetaStarD0;
+       }
+       else if (mcPart->GetPdgCode() == -421){  // D0bar{
+               AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
+               cosThetaStar = cosThetaStarD0bar;
+       }
+       else{
+               AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
+               return vectorMC;
+       }
+       if (cosThetaStar < -1 || cosThetaStar > 1) {
+               AliWarning("Invalid value for cosine Theta star, returning");
+               return isOk;
+       }
+
+       // calculate cos(Theta*) according to the method implemented herein
+
+       Double_t cts = 9999.;
+       cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
+       if (cts < -1 || cts > 1) {
+               AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVar method");
+       }
+       if (TMath::Abs(cts - cosThetaStar)>0.001) {
+               AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
+       }
+       
+       Double_t cT = decay->Ct(421);
+
+       // calculate cT from the method implemented herein
+       Double_t cT1 = 0.;
+       cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
+
+       if (TMath::Abs(cT1 - cT)>0.001) {
+               AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
+       }
+       
+       // get the pT of the daughters
+       
+       Double_t pTpi = 0.;
+       Double_t pTK = 0.;
+       
+       if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
+               pTpi = mcPartDaughter0->Pt();
+               pTK = mcPartDaughter1->Pt();
+       }
+       else {
+               pTpi = mcPartDaughter1->Pt();
+               pTK = mcPartDaughter0->Pt();
+       }
+
+       vectorMC[0] = mcPart->Pt();
+       vectorMC[1] = mcPart->Y() ;
+       vectorMC[2] = cosThetaStar ;
+       vectorMC[3] = pTpi ;
+       vectorMC[4] = pTK ;
+       vectorMC[5] = cT*1.E4 ;  // in micron
+       isOk = kTRUE;
+       return isOk;
+}
+
diff --git a/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.h b/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.h
new file mode 100644 (file)
index 0000000..9a7449d
--- /dev/null
@@ -0,0 +1,81 @@
+#ifndef ALICFHEAVYFLAVOURTASKMULTIVAR_H
+#define ALICFHEAVYFLAVOURTASKMULTIVAR_H
+/**************************************************************************
+ * 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 for HF corrections as a function of many variables
+// Author : C. Zampolli, CERN
+//-----------------------------------------------------------------------
+
+
+#include "AliAnalysisTaskSE.h"
+
+class TH1I;
+class TParticle ;
+class TFile ;
+class TClonesArray ;
+class AliCFManager;
+class AliAODRecoDecay;
+class AliAODRecoDecayHF2Prong;
+class AliAODMCParticle;
+
+class AliCFHeavyFlavourTaskMultiVar : public AliAnalysisTaskSE {
+  public:
+
+  enum {
+    kStepGenerated       = 0,
+    kStepReconstructed   = 1,
+  };
+
+  AliCFHeavyFlavourTaskMultiVar();
+  AliCFHeavyFlavourTaskMultiVar(const Char_t* name);
+  AliCFHeavyFlavourTaskMultiVar& operator= (const AliCFHeavyFlavourTaskMultiVar& c);
+  AliCFHeavyFlavourTaskMultiVar(const AliCFHeavyFlavourTaskMultiVar& c);
+  virtual ~AliCFHeavyFlavourTaskMultiVar();
+
+  // ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects
+  void     UserCreateOutputObjects();
+  void     UserExec(Option_t *option);
+  void     Terminate(Option_t *);
+  
+  // CORRECTION FRAMEWORK RELATED FUNCTIONS
+  void           SetCFManager(AliCFManager* io) {fCFManager = io;}   // global correction manager
+  AliCFManager * GetCFManager()                 {return fCFManager;} // get corr manager
+
+  void     SetPDG(Int_t code) {fPDG = code; }     // defines the PDG code of searched HF
+  Int_t    IsMcVtx(AliAODRecoDecayHF2Prong* vtx); // checks if the AliAODv0 can be associated, returns mother label
+  Int_t    GetVtxLabel(Int_t* labels);            // returns label of vertex given the daughter labels
+  Double_t CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const;  // returns cos(ThetaStar) of the D0 decay
+  Double_t CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const;            // returns cT of the D0 decay
+  void     SetFillFromGenerated(Bool_t flag) {fFillFromGenerated = flag;}
+  Bool_t   GetFillFromGenerated() const {return fFillFromGenerated;}
+  Bool_t   GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const;
+
+ protected:
+  Int_t           fPDG;         //  PDG code of searched V0's
+  AliCFManager   *fCFManager  ; //  pointer to the CF manager
+  TH1I *fHistEventsProcessed;   //! simple histo for monitoring the number of events processed
+  Int_t fCountMC;               //  MC particle found
+  Int_t fEvents;                //  n. of events
+  Bool_t fFillFromGenerated;    //  flag to indicate whether data container should be filled 
+                                //  with generated values also for reconstructed particles
+  Int_t fSkipped;               //  n. of reco decays skipped after checking on the decay (MUST be D0->K+pi)
+
+  
+  ClassDef(AliCFHeavyFlavourTaskMultiVar,0); // class for HF corrections as a function of many variables
+};
+
+#endif
diff --git a/PWG3/vertexingHF/ReadCFHeavyFlavourOutput.C b/PWG3/vertexingHF/ReadCFHeavyFlavourOutput.C
new file mode 100644 (file)
index 0000000..1c8368a
--- /dev/null
@@ -0,0 +1,766 @@
+#include <Riostream.h>
+
+extern TRandom *gRandom;
+extern TBenchmark *gBenchmark;
+extern TSystem *gSystem;
+
+void ReadCFHeavyFlavourOutput(){
+
+       // example macro for reading a container coming from a CF analysis
+       // the generated and reconstructed distributions are got, 
+       // the efficiencies are calculated,
+       // the reco distributions are corrected,
+       // the results plotted
+       
+       gROOT->SetStyle("Plain");
+       gStyle->SetPalette(1);
+       gStyle->SetOptStat(1110);
+       gStyle->SetPalette(1);
+       gStyle->SetCanvasColor(0);
+       gStyle->SetFrameFillColor(0);
+       gStyle->SetOptTitle(0);
+       
+       gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include  -I$ROOTSYS/include");
+       gSystem->Load("libANALYSIS.so");
+       gSystem->Load("libANALYSISalice.so");
+       gSystem->Load("$ALICE_ROOT/CORRFW/libCORRFW.so") ;
+       
+       //Setting up the container grid... 
+       
+       const Int_t nstep=2; //number of selection steps  (just 2 in this ex)
+       
+       const Int_t nvar=6; //number of variables on the grid: pt, y, cosThetaStar, ptPi, ptK, cT
+       
+       // Flag the sel steps. In this example, we have two, may be any nstep
+       Int_t stepGen=0;
+       Int_t stepRec=1;
+       
+       // the sensitive variables, their indeces
+       UInt_t ipt  =0;
+       UInt_t iy   =1;
+       UInt_t icTS = 2; // cTS stands for cosThetaStar
+       UInt_t iptPi= 3;
+       UInt_t iptK = 4;
+       UInt_t icT  = 5;
+       
+       // Read the  container from file
+       TFile *file = new TFile("output.root");
+       AliCFContainer *data = (AliCFContainer*) (file->Get("container"));
+       
+       // Make some 1 & 2-D projections..
+       // MC level
+       TCanvas *cmc1 = new TCanvas("cmc1","The MC distributions",0,800,900,1200);
+       TCanvas *cmc2 = new TCanvas("cmc2","The MC distributions",0,800,900,1200);
+       TCanvas *cmcpt = new TCanvas("cmcpt","pt distribution from MC",50,50,550,550);
+       TCanvas *cmcy = new TCanvas("cmcy","y distribution from MC",50,50,550,550);
+       TCanvas *cmccTS = new TCanvas("cmcTS","cosThetaStar distribution from MC",50,50,550,550);  // cTS stands for cosThetaStar
+       TCanvas *cmcptPi = new TCanvas("cmcptPi","pt_pi distribution from MC",50,50,550,550);
+       TCanvas *cmcptK = new TCanvas("cmcptK","pt_K distribution from MC",50,50,550,550);
+       TCanvas *cmccT = new TCanvas("cmccT","cT distribution from MC",50,50,550,550);
+
+       // Reco-aod level
+       TCanvas *cpt = new TCanvas("cpt","pt distribution from reco aod",50,50,550,550);
+       TCanvas *cy = new TCanvas("cy","y distribution from reco aod",50,50,550,550);
+       TCanvas *ccTS = new TCanvas("cTS","cosThetaStar distribution from reco aod",50,50,550,550); 
+       TCanvas *cptPi = new TCanvas("cptPi","pt_pi distribution from reco aod",50,50,550,550);
+       TCanvas *cptK = new TCanvas("cptK","pt_K distribution from reco aod",50,50,550,550);
+       TCanvas *ccT = new TCanvas("ccT","cT distribution from reco aod",50,50,550,550);
+
+       // some 2D distributions in pt and y
+       TCanvas *cGen2d = new TCanvas("cGen2d","2D distribution from MC ",50,50,550,550);
+       TCanvas *cRec2d = new TCanvas("cRec2d","2D distribution from reco aod",50,50,550,550);
+       
+       // 2D plots 
+       TH2D *hMCpty2d = data->ShowProjection(ipt, iy, stepGen);
+       hMCpty2d->SetLineColor(2);
+       hMCpty2d->SetLineWidth(3);
+       hMCpty2d->SetMarkerColor(2);
+       hMCpty2d->SetMarkerStyle(20);
+       hMCpty2d->GetXaxis()->SetTitleOffset(1.2);
+       hMCpty2d->GetYaxis()->SetTitleOffset(1.5);
+       hMCpty2d->GetXaxis()->SetTitle("p_{T} (GeV/c), MC data");
+       hMCpty2d->GetYaxis()->SetTitle("y, MC data");
+       hMCpty2d->Draw("text");
+       cGen2d->cd();
+       cGen2d->SetLeftMargin(0.15);
+       cGen2d->SetRightMargin(0.05);
+       cGen2d->Update();  
+
+       TH2D *hRECpty2d = data->ShowProjection(ipt, iy, stepRec);
+       hRECpty2d->Sumw2();
+       hRECpty2d->SetLineColor(4);
+       hRECpty2d->SetLineWidth(3);
+       hRECpty2d->SetMarkerColor(4);
+       hRECpty2d->SetMarkerStyle(20);
+       hRECpty2d->GetXaxis()->SetTitleOffset(1.2);
+       hRECpty2d->GetYaxis()->SetTitleOffset(1.5);
+       hRECpty2d->GetXaxis()->SetTitle("p_{T} (GeV/c), AOD");
+       hRECpty2d->GetYaxis()->SetTitle("y, AOD");
+       hRECpty2d->Draw("text");
+       cRec2d->cd();
+       cRec2d->SetLeftMargin(0.15);
+       cRec2d->SetRightMargin(0.05);
+       cRec2d->Update();  
+
+       // MC + REC 1D plots 
+       // pt, y, cosThetaStar
+       cmc1->Divide(2,3);   
+
+       cmc1->cd(1);
+       TH1D *hMCpt1D = data->ShowProjection(ipt, stepGen);
+       Double_t maxpt = hMCpt1D->GetMaximum();
+       hMCpt1D->GetYaxis()->SetRangeUser(0,maxpt*1.2);
+       hMCpt1D->Sumw2();
+       hMCpt1D->SetMinimum(0.01);
+       hMCpt1D->SetLineColor(2);
+       hMCpt1D->SetLineWidth(3);
+       hMCpt1D->SetMarkerColor(2);
+       hMCpt1D->SetFillColor(2);
+       hMCpt1D->SetFillStyle(3005);
+       hMCpt1D->SetMarkerStyle(20);
+       hMCpt1D->GetXaxis()->SetTitleOffset(1.2);
+       hMCpt1D->GetXaxis()->SetTitle("p_{T} (GeV/c), MC data");
+       hMCpt1D->Draw("hist");
+       cmcpt->cd();
+       cmcpt->SetLeftMargin(0.15);
+       cmcpt->SetRightMargin(0.05);
+       hMCpt1D->Draw("hist");
+       hMCpt1D->Draw("err same");
+       cmcpt->Update();
+
+       cmc1->cd(2);
+       TH1D *hRECpt1D = data->ShowProjection(ipt, stepRec);
+       hRECpt1D->GetYaxis()->SetRangeUser(0,maxpt*1.2);
+       hRECpt1D->SetLineColor(4);
+       hRECpt1D->SetLineWidth(3);
+       hRECpt1D->SetMarkerColor(4);
+       hRECpt1D->SetFillColor(4);
+       hRECpt1D->SetFillStyle(3004);
+       hRECpt1D->SetMarkerStyle(20);
+       hRECpt1D->GetXaxis()->SetTitleOffset(1.2);
+       hRECpt1D->Sumw2();
+       hRECpt1D->SetMinimum(0.01);
+       hRECpt1D->GetXaxis()->SetTitle("p_{T} (GeV/c), AOD");
+       hRECpt1D->Draw("hist");
+       cpt->cd();
+       cpt->SetLeftMargin(0.15);
+       cpt->SetRightMargin(0.05);
+       hRECpt1D->Draw("hist");
+       hRECpt1D->Draw("err same");
+       cpt->Update();
+
+       cmc1->cd(3);
+       TH1D *hMCy1D = data->ShowProjection(iy, stepGen);
+       Double_t maxy = hMCy1D->GetMaximum();
+       hMCy1D->GetYaxis()->SetRangeUser(0,maxy*1.2);
+       hMCy1D->SetLineColor(2);
+       hMCy1D->SetLineWidth(3);
+       hMCy1D->SetMarkerColor(2);
+       hMCy1D->SetFillColor(2);
+       hMCy1D->SetFillStyle(3005);
+       hMCy1D->SetMarkerStyle(20);
+       hMCy1D->GetXaxis()->SetTitleOffset(1.2);
+       hMCy1D->Sumw2();
+       hMCy1D->SetMinimum(0.01);
+       hMCy1D->GetXaxis()->SetTitle("y, MC data");
+       hMCy1D->Draw("hist");
+       cmcy->cd();
+       cmcy->SetLeftMargin(0.15);
+       cmcy->SetRightMargin(0.05);
+       hMCy1D->Draw("hist");
+       hMCy1D->Draw("err same");
+       cmcy->Update();
+
+       cmc1->cd(4);
+       TH1D *hRECy1D = data->ShowProjection(iy, stepRec);
+       hRECy1D->GetYaxis()->SetRangeUser(0,maxy*1.2);
+       hRECy1D->SetLineColor(4);
+       hRECy1D->SetLineWidth(3);
+       hRECy1D->SetMarkerColor(4);
+       hRECy1D->SetFillColor(4);
+       hRECy1D->SetFillStyle(3004);
+       hRECy1D->SetMarkerStyle(20);
+       hRECy1D->GetXaxis()->SetTitleOffset(1.2);
+       hRECy1D->Sumw2();
+       hRECy1D->SetMinimum(0.01);
+       hRECy1D->Draw("hist");
+       cy->cd();
+       cy->SetLeftMargin(0.15);
+       cy->SetRightMargin(0.05);
+       hRECy1D->GetXaxis()->SetTitle("y, AOD");
+       hRECy1D->Draw("hist");
+       hRECy1D->Draw("err same");
+       cy->Update();
+
+       cmc1->cd(5);
+       TH1D *hMCcTS1D = data->ShowProjection(icTS, stepGen);
+       Double_t maxcTS = hMCcTS1D->GetMaximum();
+       hMCcTS1D->GetYaxis()->SetRangeUser(0,maxcTS*1.2);
+       hMCcTS1D->SetLineColor(2);
+       hMCcTS1D->SetLineWidth(3);
+       hMCcTS1D->SetMarkerColor(2);
+       hMCcTS1D->SetFillColor(2);
+       hMCcTS1D->SetFillStyle(3005);
+       hMCcTS1D->SetMarkerStyle(20);
+       hMCcTS1D->GetXaxis()->SetTitleOffset(1.2);
+       hMCcTS1D->Sumw2();
+       hMCcTS1D->SetMinimum(0.01);
+       hMCcTS1D->GetXaxis()->SetTitle("cosThetaStar, MC data");
+       hMCcTS1D->Draw("hist");
+       cmccTS->cd();
+       cmccTS->SetLeftMargin(0.15);
+       cmccTS->SetRightMargin(0.05);
+       hMCcTS1D->Draw("hist");
+       hMCcTS1D->Draw("err same");
+       cmccTS->Update();
+
+       cmc1->cd(6);
+       TH1D *hRECcTS1D = data->ShowProjection(icTS, stepRec);
+       hRECcTS1D->GetYaxis()->SetRangeUser(0,maxcTS*1.2);
+       hRECcTS1D->SetLineColor(4);
+       hRECcTS1D->SetLineWidth(3);
+       hRECcTS1D->SetMarkerColor(4);
+       hRECcTS1D->SetFillColor(4);
+       hRECcTS1D->SetFillStyle(3004);
+       hRECcTS1D->SetMarkerStyle(20);
+       hRECcTS1D->GetXaxis()->SetTitleOffset(1.2);
+       hRECcTS1D->Sumw2();
+       hRECcTS1D->SetMinimum(0.01);
+       hRECcTS1D->Draw("hist");
+       ccTS->cd();
+       ccTS->SetLeftMargin(0.15);
+       ccTS->SetRightMargin(0.05);
+       hRECcTS1D->GetXaxis()->SetTitle("cosThetaStar, AOD");
+       hRECcTS1D->Draw("hist");
+       hRECcTS1D->Draw("err same");
+       ccTS->Update();
+
+       // ptPi, ptK, cT
+       cmc2->Divide(2,3); 
+
+       cmc2->cd(1);
+       TH1D *hMCptPi1D = data->ShowProjection(iptPi, stepGen);
+       Double_t maxptPi = hMCptPi1D->GetMaximum();
+       hMCptPi1D->GetYaxis()->SetRangeUser(0,maxptPi*1.2);
+       hMCptPi1D->Sumw2();
+       hMCptPi1D->SetMinimum(0.01);
+       hMCptPi1D->SetLineColor(2);
+       hMCptPi1D->SetLineWidth(3);
+       hMCptPi1D->SetMarkerColor(2);
+       hMCptPi1D->SetFillColor(2);
+       hMCptPi1D->SetFillStyle(3005);
+       hMCptPi1D->SetMarkerStyle(20);
+       hMCptPi1D->GetXaxis()->SetTitleOffset(1.2);
+       hMCptPi1D->GetXaxis()->SetTitle("p_{T, #pi} (GeV/c), MC data");
+       hMCptPi1D->Draw("hist");
+       cmcptPi->cd();
+       cmcptPi->SetLeftMargin(0.15);
+       cmcptPi->SetRightMargin(0.05);
+       hMCptPi1D->Draw("hist");
+       hMCptPi1D->Draw("err same");
+       cmcptPi->Update();
+
+       cmc2->cd(2);
+       TH1D *hRECptPi1D = data->ShowProjection(iptPi, stepRec);
+       hRECptPi1D->GetYaxis()->SetRangeUser(0,maxptPi*1.2);
+       hRECptPi1D->SetLineColor(4);
+       hRECptPi1D->SetLineWidth(3);
+       hRECptPi1D->SetMarkerColor(4);
+       hRECptPi1D->SetFillColor(4);
+       hRECptPi1D->SetFillStyle(3004);
+       hRECptPi1D->SetMarkerStyle(20);
+       hRECptPi1D->GetXaxis()->SetTitleOffset(1.2);
+       hRECptPi1D->Sumw2();
+       hRECptPi1D->SetMinimum(0.01);
+       hRECptPi1D->GetXaxis()->SetTitle("p_{T, #pi} (GeV/c), AOD");
+       hRECptPi1D->Draw("hist");
+       cptPi->cd();
+       cptPi->SetLeftMargin(0.15);
+       cptPi->SetRightMargin(0.05);
+       hRECptPi1D->Draw("hist");
+       hRECptPi1D->Draw("err same");
+       cptPi->Update();
+
+       cmc2->cd(3);
+       TH1D *hMCptK1D = data->ShowProjection(iptK, stepGen);
+       Double_t maxptK = hMCptK1D->GetMaximum();
+       hMCptK1D->GetYaxis()->SetRangeUser(0,maxptK*1.2);
+       hMCptK1D->SetLineColor(2);
+       hMCptK1D->SetLineWidth(3);
+       hMCptK1D->SetMarkerColor(2);
+       hMCptK1D->SetFillColor(2);
+       hMCptK1D->SetFillStyle(3005);
+       hMCptK1D->SetMarkerStyle(20);
+       hMCptK1D->GetXaxis()->SetTitleOffset(1.2);
+       hMCptK1D->Sumw2();
+       hMCptK1D->SetMinimum(0.01);
+       hMCptK1D->GetXaxis()->SetTitle("p_{T, K} (GeV/c), MC data");
+       hMCptK1D->Draw("hist");
+       cmcptK->cd();
+       cmcptK->SetLeftMargin(0.15);
+       cmcptK->SetRightMargin(0.05);
+       hMCptK1D->Draw("hist");
+       hMCptK1D->Draw("err same");
+       cmcptK->Update();
+
+       cmc2->cd(4);
+       TH1D *hRECptK1D = data->ShowProjection(iptK, stepRec);
+       hRECptK1D->GetYaxis()->SetRangeUser(0,maxptK*1.2);
+       hRECptK1D->SetLineColor(4);
+       hRECptK1D->SetLineWidth(3);
+       hRECptK1D->SetMarkerColor(4);
+       hRECptK1D->SetFillColor(4);
+       hRECptK1D->SetFillStyle(3004);
+       hRECptK1D->SetMarkerStyle(20);
+       hRECptK1D->GetXaxis()->SetTitleOffset(1.2);
+       hRECptK1D->Sumw2();
+       hRECptK1D->SetMinimum(0.01);
+       hRECptK1D->Draw("hist");
+       cptK->cd();
+       cptK->SetLeftMargin(0.15);
+       cptK->SetRightMargin(0.05);
+       hRECptK1D->GetXaxis()->SetTitle("p_{T, K} (GeV/c), AOD");
+       hRECptK1D->Draw("hist");
+       hRECptK1D->Draw("err same");
+       cptK->Update();
+
+       cmc2->cd(5);
+       TH1D *hMCcT1D = data->ShowProjection(icT, stepGen);
+       Double_t maxcT = hMCcT1D->GetMaximum();
+       hMCcT1D->GetYaxis()->SetRangeUser(0,maxcT*1.2);
+       hMCcT1D->SetLineColor(2);
+       hMCcT1D->SetLineWidth(3);
+       hMCcT1D->SetMarkerColor(2);
+       hMCcT1D->SetFillColor(2);
+       hMCcT1D->SetFillStyle(3005);
+       hMCcT1D->SetMarkerStyle(20);
+       hMCcT1D->GetXaxis()->SetTitleOffset(1.2);
+       hMCcT1D->Sumw2();
+       hMCcT1D->SetMinimum(0.01);
+       hMCcT1D->GetXaxis()->SetTitle("ct (#mum), MC data");
+       hMCcT1D->Draw("hist");
+       cmccT->cd();
+       cmccT->SetLeftMargin(0.15);
+       cmccT->SetRightMargin(0.05);
+       hMCcT1D->Draw("hist");
+       hMCcT1D->Draw("err same");
+       cmccT->Update();
+
+       cmc2->cd(6);
+       TH1D *hRECcT1D = data->ShowProjection(icT, stepRec);
+       hRECcT1D->GetYaxis()->SetRangeUser(0,maxcT*1.2);
+       hRECcT1D->SetLineColor(4);
+       hRECcT1D->SetLineWidth(3);
+       hRECcT1D->SetMarkerColor(4);
+       hRECcT1D->SetFillColor(4);
+       hRECcT1D->SetFillStyle(3004);
+       hRECcT1D->SetMarkerStyle(20);
+       hRECcT1D->GetXaxis()->SetTitleOffset(1.2);
+       hRECcT1D->Sumw2();
+       hRECcT1D->SetMinimum(0.01);
+       hRECcT1D->Draw("hist");
+       ccT->cd();
+       ccT->SetLeftMargin(0.15);
+       ccT->SetRightMargin(0.05);
+       hRECcT1D->GetXaxis()->SetTitle("c#t (#mum), AOD");
+       hRECcT1D->Draw("hist");
+       hRECcT1D->Draw("err same");
+       ccT->Update();
+
+       /*
+       // printing on eps files
+       cmc1->Print("Plots/dataMC_pt_y_cTS.gif");
+       cmc2->Print("Plots/dataMC_ptPi_ptK_cT.gif");
+       cmcpt->Print("Plots/pt_Gen.eps");
+       cmcy->Print("Plots/y_Gen.eps");
+       cmccTS->Print("Plots/cTS_Gen.eps");
+       cmcptPi->Print("Plots/ptPi_Gen.eps");
+       cmcptK->Print("Plots/ptK_Gen.eps");
+       cmccT->Print("Plots/cT_Gen.eps");
+       cpt->Print("Plots/pt_Rec.eps");
+       cy->Print("Plots/y_Rec.eps");
+       ccTS->Print("Plots/cTS_Rec.eps");
+       cptPi->Print("Plots/ptPi_Rec.eps");
+       cptK->Print("Plots/ptK_Rec.eps");
+       ccT->Print("Plots/cT_Rec.eps");
+       cGen2d->Print("Plots/pt_y_Gen_2D.eps");  
+       cRec2d->Print("Plots/pt_y_Rec_2D.eps");  
+
+       // printing on gif files
+       cmc1->Print("Plots/dataMC_pt_y_cTS.gif");
+       cmc2->Print("Plots/dataMC_ptPi_ptK_cT.gif");
+       cmcpt->Print("Plots/pt_Gen.eps");
+       cmcy->Print("Plots/y_Gen.eps");
+       cmccTS->Print("Plots/cTS_Gen.eps");
+       cmcptPi->Print("Plots/ptPi_Gen.eps");
+       cmcptK->Print("Plots/ptK_Gen.eps");
+       cmccT->Print("Plots/cT_Gen.eps");
+       cpt->Print("Plots/pt_Rec.eps");
+       cy->Print("Plots/y_Rec.eps");
+       ccTS->Print("Plots/cTS_Rec.eps");
+       cptPi->Print("Plots/ptPi_Rec.eps");
+       cptK->Print("Plots/ptK_Rec.eps");
+       ccT->Print("Plots/cT_Rec.eps");
+       cGen2d->Print("Plots/pt_y_Gen_2D.eps");  
+       cRec2d->Print("Plots/pt_y_Rec_2D.eps");  
+       */
+
+       //construct the efficiency grid from the data container 
+       AliCFEffGrid *eff = new AliCFEffGrid("eff"," The efficiency",*data);
+       eff->CalculateEfficiency(stepRec,stepGen); //eff= step1/step0
+
+       //The efficiency along the variables, and some 2-D projections
+       TCanvas *ceff =new TCanvas("ceff"," Efficiency",0,0,1600,1200);
+       ceff->Divide(3,2);
+       TCanvas *ceff2D =new TCanvas("ceff2D"," Efficiency for pt and y",50,50,550,550);
+       TCanvas *ceffpt = new TCanvas("ceffpt","Efficiency vs pt",50,50,550,550);
+       TCanvas *ceffy = new TCanvas("ceffy","Efficiency vs y",50,50,550,550);
+       TCanvas *ceffcTS = new TCanvas("ceffcTS","Efficiency vs cosThetaStar",50,50,550,550);
+       TCanvas *ceffptPi = new TCanvas("ceffptPi","Efficiency vs ptPi",50,50,550,550);
+       TCanvas *ceffptK = new TCanvas("ceffptK","Efficiency vs ptK",50,50,550,550);
+       TCanvas *ceffcT = new TCanvas("ceffcT","Efficiency vs cT",50,50,550,550);
+       TCanvas *ceff2Dtext = new TCanvas("ceff2Dtext","Text plot for efficiency in pt and y",50,50,550,550);
+
+       ceff->cd(1);
+       TH1D *hpteffCF = eff->Project(ipt); //the efficiency vs pt
+       hpteffCF->Sumw2();
+       //hpteffCF->SetMinimum(0.01);
+       hpteffCF->SetLineColor(8);
+       hpteffCF->SetLineWidth(3);
+       hpteffCF->SetMarkerColor(8);
+       hpteffCF->SetMarkerStyle(20);
+       hpteffCF->GetXaxis()->SetTitleOffset(1.2);
+       hpteffCF->GetYaxis()->SetTitleOffset(1.5);
+       hpteffCF->Draw("hist");
+       ceffpt->cd();
+       ceffpt->SetLeftMargin(0.15);
+       ceffpt->SetRightMargin(0.05);
+       hpteffCF->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+       hpteffCF->GetYaxis()->SetTitle("Efficiency");
+       hpteffCF->Draw("hist");
+       hpteffCF->Draw("err same");
+
+       ceff->cd(2);
+       TH1D *hyeffCF = eff->Project(iy); //the efficiency vs y
+       hyeffCF->Sumw2();
+       //hyeffCF->SetMinimum(0.01);
+       hyeffCF->SetLineColor(8);
+       hyeffCF->SetLineWidth(3);
+       hyeffCF->SetMarkerColor(8);
+       hyeffCF->SetMarkerStyle(20);
+       hyeffCF->GetXaxis()->SetTitleOffset(1.2);
+       hyeffCF->GetYaxis()->SetTitleOffset(1.5);
+       hyeffCF->Draw("hist");
+       ceffy->cd();
+       ceffy->SetLeftMargin(0.15);
+       ceffy->SetRightMargin(0.05);
+       hyeffCF->GetXaxis()->SetTitle("y");
+       hyeffCF->GetYaxis()->SetTitle("Efficiency");
+       hyeffCF->Draw("hist");
+       hyeffCF->Draw("err same");
+
+       ceff->cd(3);
+       TH1D *hcTSeffCF = eff->Project(icTS); //the efficiency vs cosThetaStar
+       hcTSeffCF->Sumw2();
+       //hcTSeffCF->SetMinimum(0.01);
+       hcTSeffCF->SetLineColor(8);
+       hcTSeffCF->SetLineWidth(3);
+       hcTSeffCF->SetMarkerColor(8);
+       hcTSeffCF->SetMarkerStyle(20);
+       hcTSeffCF->GetXaxis()->SetTitleOffset(1.2);
+       hcTSeffCF->GetYaxis()->SetTitleOffset(1.5);
+       hcTSeffCF->Draw("hist");
+       ceffcTS->cd();
+       ceffcTS->SetLeftMargin(0.15);
+       ceffcTS->SetRightMargin(0.05);
+       hcTSeffCF->GetXaxis()->SetTitle("cosThetaStar");
+       hcTSeffCF->GetYaxis()->SetTitle("Efficiency");
+       hcTSeffCF->Draw("hist");
+       hcTSeffCF->Draw("err same");
+
+       ceff->cd(4);
+       TH1D *hptPieffCF = eff->Project(iptPi); //the efficiency vs ptPi
+       hptPieffCF->Sumw2();
+       //hptPieffCF->SetMinimum(0.01);
+       hptPieffCF->SetLineColor(8);
+       hptPieffCF->SetLineWidth(3);
+       hptPieffCF->SetMarkerColor(8);
+       hptPieffCF->SetMarkerStyle(20);
+       hptPieffCF->GetXaxis()->SetTitleOffset(1.2);
+       hptPieffCF->GetYaxis()->SetTitleOffset(1.5);
+       hptPieffCF->Draw("hist");
+       ceffptPi->cd();
+       ceffptPi->SetLeftMargin(0.15);
+       ceffptPi->SetRightMargin(0.05);
+       hptPieffCF->GetXaxis()->SetTitle("p_{T, #pi} (GeV/c)");
+       hptPieffCF->GetYaxis()->SetTitle("Efficiency");
+       hptPieffCF->Draw("hist");
+       hptPieffCF->Draw("err same");
+
+       ceff->cd(5);
+       TH1D *hptKeffCF = eff->Project(iptK); //the efficiency vs ptK
+       hptKeffCF->Sumw2();
+       //hptKeffCF->SetMinimum(0.01);
+       hptKeffCF->SetLineColor(8);
+       hptKeffCF->SetLineWidth(3);
+       hptKeffCF->SetMarkerColor(8);
+       hptKeffCF->SetMarkerStyle(20);
+       hptKeffCF->GetXaxis()->SetTitleOffset(1.2);
+       hptKeffCF->GetYaxis()->SetTitleOffset(1.5);
+       hptKeffCF->Draw("hist");
+       ceffptK->cd();
+       ceffptK->SetLeftMargin(0.15);
+       ceffptK->SetRightMargin(0.05);
+       hptKeffCF->GetXaxis()->SetTitle("p_{T, K} (GeV/c)");
+       hptKeffCF->GetYaxis()->SetTitle("Efficiency");
+       hptKeffCF->Draw("hist");
+       hptKeffCF->Draw("err same");
+
+       ceff->cd(6);
+       TH1D *hcTeffCF = eff->Project(icT); //the efficiency vs cT
+       hcTeffCF->Sumw2();
+       //hcTeffCF->SetMinimum(0.01);
+       hcTeffCF->SetLineColor(8);
+       hcTeffCF->SetLineWidth(3);
+       hcTeffCF->SetMarkerColor(8);
+       hcTeffCF->SetMarkerStyle(20);
+       hcTeffCF->GetXaxis()->SetTitleOffset(1.2);
+       hcTeffCF->GetYaxis()->SetTitleOffset(1.5);
+       hcTeffCF->Draw("hist");
+       ceffcT->cd();
+       ceffcT->SetLeftMargin(0.15);
+       ceffcT->SetRightMargin(0.05);
+       hcTeffCF->GetXaxis()->SetTitle("c#t (#mum)");
+       hcTeffCF->GetYaxis()->SetTitle("Efficiency");
+       hcTeffCF->Draw("hist");
+       hcTeffCF->Draw("err same");
+
+       ceff2D->cd();
+       TH2D *hptyeffCF = eff->Project(ipt,iy); //look at the numerator
+       //hptyeffCF->SetMinimum(0.01);
+       hptyeffCF->SetMarkerColor(8);
+       hptyeffCF->SetLineColor(8);
+       hptyeffCF->SetMinimum(0.01);
+       hptyeffCF->Draw("lego");
+       ceff2Dtext->cd();
+       hptyeffCF->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+       hptyeffCF->GetXaxis()->SetTitleOffset(1.2);
+       hptyeffCF->GetYaxis()->SetTitle("y");
+       hptyeffCF->GetYaxis()->SetTitleOffset(1.2);
+       hptyeffCF->Draw("text");
+
+       /*
+       // printing eps files
+       ceff->Print("Plots/efficiencies.eps");
+       ceffpt->Print("Plots/effpt.eps");
+       ceffy->Print("Plots/effy.eps");
+       ceffcTS->Print("Plots/effcTS.eps");
+       ceffptPi->Print("Plots/effptPi.eps");
+       ceffptK->Print("Plots/effptK.eps");
+       ceffcT->Print("Plots/effcT.eps");
+       ceff2D->Print("Plots/eff2D_pt_y.eps");
+       ceff2Dtext->Print("Plots/eff2d_text_pt_y.eps");
+
+       // printing gif files
+       ceff->Print("Plots/efficiencies.gif");
+       ceffpt->Print("Plots/effpt.gif");
+       ceffy->Print("Plots/effy.gif");
+       ceffcTS->Print("Plots/effcTS.gif");
+       ceffptPi->Print("Plots/effptPi.gif");
+       ceffptK->Print("Plots/effptK.gif");
+       ceffcT->Print("Plots/effcT.gif");
+       ceff2D->Print("Plots/eff2D_pt_y.gif");
+       ceff2Dtext->Print("Plots/eff2d_text_pt_y.gif");
+       */
+
+       // applying efficiencies - using projections, not the ApplyEffCorrection method
+
+       TCanvas *cmultpiplypt = new TCanvas("cmultiplypt","Reco From Eff in pt distribution",50,50,550,550);
+       TCanvas *cmultpiplyy = new TCanvas("cmultiplyy","Reco From Eff in y distribution",50,50,550,550);
+       TCanvas *cmultpiplycTS = new TCanvas("cmultiplycTS","Reco From Eff in cTS distribution",50,50,550,550);
+       TCanvas *cmultpiplyptPi = new TCanvas("cmultiplyptPi","Reco From Eff in ptPi distribution",50,50,550,550);
+       TCanvas *cmultpiplyptK = new TCanvas("cmultiplyptK","Reco From Eff in ptK distribution",50,50,550,550);
+       TCanvas *cmultpiplycT = new TCanvas("cmultiplycT","Reco From Eff in cT distribution",50,50,550,550);
+
+       TH1D *hmultiplypt = new TH1D("hmultiplypt","hmultiplypt",13,0,10);
+       cout << " bin for histo MC = " << hMCpt1D->GetNbinsX() << " while for RECO histo = " << hRECpt1D->GetNbinsX() << " while for efficiency histo = " << hpteffCF->GetNbinsX() << endl;
+
+       const Double_t ptmin_0_4 =  0.0 ;
+       const Double_t ptmax_0_4 =  4.0 ;
+       const Double_t ptmin_4_8 =  4.0 ;
+       const Double_t ptmax_4_8 =  8.0 ;
+       const Double_t ptmin_8_10 =  8.0 ;
+       const Double_t ptmax_8_10 =  10.0 ;
+       const Int_t nbin0_0_4  = 8 ; //bins in pt
+       const Int_t nbin0_4_8  = 4 ; //bins in pt
+       const Int_t nbin0_8_10  = 1 ; //bins in pt
+       const Int_t nbins = nbin0_0_4 + nbin0_4_8 + nbin0_8_10;
+       Double_t binLim0[nbins+1];
+       for(Int_t i=0; i<=nbin0_0_4; i++) binLim0[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin0_0_4*(Double_t)i ; 
+       if (binLim0[nbin0_0_4] != ptmin_4_8)  {
+               printf("Calculated bin lim for pt - 1st range - differs from expected!\n");
+       }
+       for(Int_t i=0; i<=nbin0_4_8; i++) binLim0[i+nbin0_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin0_4_8*(Double_t)i ; 
+       if (binLim0[nbin0_0_4+nbin0_4_8] != ptmin_8_10)  {
+               printf("Calculated bin lim for pt - 2nd range - differs from expected!\n");
+       }
+       for(Int_t i=0; i<=nbin0_8_10; i++) binLim0[i+nbin0_0_4+nbin0_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin0_8_10*(Double_t)i ; 
+
+       hmultiplypt->SetBins(nbins,binLim0);
+       hmultiplypt->Divide(hRECpt1D,hpteffCF,1,1,"b");
+       hmultiplypt->SetMinimum(0);
+       hmultiplypt->SetMaximum(hMCpt1D->GetMaximum());
+       hmultiplypt->SetLineColor(38);
+       hmultiplypt->SetLineWidth(3);
+       hmultiplypt->SetMarkerColor(38);
+       hmultiplypt->SetMarkerStyle(20);
+       hmultiplypt->GetXaxis()->SetTitleOffset(1.2);
+       hmultiplypt->GetYaxis()->SetTitleOffset(1.5);
+       hmultiplypt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+       cmultiplypt->SetLeftMargin(0.15);
+       cmultiplypt->SetRightMargin(0.05);
+       hmultiplypt->SetFillStyle(3004);
+       hmultiplypt->SetFillColor(38);
+       cmultiplypt->cd();
+       hmultiplypt->Draw("hist err");
+       hMCpt1D->Draw("hist same err");
+
+       TH1D *hmultiplyy = new TH1D("hmultiplyy","hmultiplyy",8,-2,2);
+       hmultiplyy->Divide(hRECy1D,hyeffCF,1,1,"b");
+       hmultiplyy->SetMinimum(0);
+       hmultiplyy->SetMaximum(hMCy1D->GetMaximum());
+       hmultiplyy->SetLineColor(38);
+       hmultiplyy->SetLineWidth(3);
+       hmultiplyy->SetMarkerColor(38);
+       hmultiplyy->SetMarkerStyle(20);
+       hmultiplyy->GetXaxis()->SetTitleOffset(1.2);
+       hmultiplyy->GetXaxis()->SetTitle("y");
+       hmultiplyy->GetYaxis()->SetTitleOffset(1.5);
+       cmultiplyy->SetLeftMargin(0.15);
+       cmultiplyy->SetRightMargin(0.05);
+       cmultiplyy->cd();
+       hmultiplyy->SetFillStyle(3004);
+       hmultiplyy->SetFillColor(38);
+       hmultiplyy->Draw("hist err");
+       hMCy1D->Draw("hist same err");
+
+       TH1D *hmultiplycTS = new TH1D("hmultiplycTS","hmultiplycTS",8,-1,1);
+       hmultiplycTS->Divide(hRECcTS1D,hcTSeffCF,1,1,"b");
+       hmultiplycTS->SetMinimum(0);
+       hmultiplycTS->SetMaximum(hMCcTS1D->GetMaximum());
+       hmultiplycTS->SetLineColor(38);
+       hmultiplycTS->SetLineWidth(3);
+       hmultiplycTS->SetMarkerColor(38);
+       hmultiplycTS->SetMarkerStyle(20);
+       hmultiplycTS->GetXaxis()->SetTitleOffset(1.2);
+       hmultiplycTS->GetXaxis()->SetTitle("cosThetaStar");
+       hmultiplycTS->GetYaxis()->SetTitleOffset(1.5);
+       cmultiplycTS->SetLeftMargin(0.15);
+       cmultiplycTS->SetRightMargin(0.05);
+       cmultiplycTS->cd();
+       hmultiplycTS->SetFillStyle(3004);
+       hmultiplycTS->SetFillColor(38);
+       hmultiplycTS->Draw("hist err");
+       hMCcTS1D->Draw("hist same err");
+
+       TH1D *hmultiplyptPi = new TH1D("hmultiplyptPi","hmultiplyptPi",13,0,10);
+       hmultiplyptPi->SetBins(nbins,binLim0);
+       hmultiplyptPi->Divide(hRECptPi1D,hptPieffCF,1,1,"b");
+       hmultiplyptPi->SetMinimum(0);
+       hmultiplyptPi->SetMaximum(hMCptPi1D->GetMaximum());
+       hmultiplyptPi->SetLineColor(38);
+       hmultiplyptPi->SetLineWidth(3);
+       hmultiplyptPi->SetMarkerColor(38);
+       hmultiplyptPi->SetMarkerStyle(20);
+       hmultiplyptPi->GetXaxis()->SetTitleOffset(1.2);
+       hmultiplyptPi->GetXaxis()->SetTitle("p_{T, #pi} (GeV/c)");
+       hmultiplyptPi->GetYaxis()->SetTitleOffset(1.5);
+       cmultiplyptPi->SetLeftMargin(0.15);
+       cmultiplyptPi->SetRightMargin(0.05);
+       cmultiplyptPi->cd();
+       hmultiplyptPi->SetFillStyle(3004);
+       hmultiplyptPi->SetFillColor(38);
+       hmultiplyptPi->Draw("hist err");
+       hMCptPi1D->Draw("hist same err");
+
+       TH1D *hmultiplyptK = new TH1D("hmultiplyptK","hmultiplyptK",13,0,10);
+       hmultiplyptK->SetBins(nbins,binLim0);
+       hmultiplyptK->Divide(hRECptK1D,hptKeffCF,1,1,"b");
+       hmultiplyptK->SetMinimum(0);
+       hmultiplyptK->SetMaximum(hMCptK1D->GetMaximum());
+       hmultiplyptK->SetLineColor(38);
+       hmultiplyptK->SetLineWidth(3);
+       hmultiplyptK->SetMarkerColor(38);
+       hmultiplyptK->SetMarkerStyle(20);
+       hmultiplyptK->GetXaxis()->SetTitleOffset(1.2);
+       hmultiplyptK->GetXaxis()->SetTitle("p_{T, K} (GeV/c)");
+       hmultiplyptK->GetYaxis()->SetTitleOffset(1.5);
+       cmultiplyptK->SetLeftMargin(0.15);
+       cmultiplyptK->SetRightMargin(0.05);
+       cmultiplyptK->cd();
+       hmultiplyptK->SetFillStyle(3004);
+       hmultiplyptK->SetFillColor(38);
+       hmultiplyptK->Draw("hist err");
+       hMCptK1D->Draw("hist same err");
+
+       TH1D *hmultiplycT = new TH1D("hmultiplycT","hmultiplycT",24,0,500);
+       hmultiplycT->Divide(hRECcT1D,hcTeffCF,1,1,"b");
+       hmultiplycT->SetMinimum(0);
+       hmultiplycT->SetMaximum(hMCcT1D->GetMaximum());
+       hmultiplycT->SetLineColor(38);
+       hmultiplycT->SetLineWidth(3);
+       hmultiplycT->SetMarkerColor(38);
+       hmultiplycT->SetMarkerStyle(20);
+       hmultiplycT->GetXaxis()->SetTitleOffset(1.2);
+       hmultiplycT->GetXaxis()->SetTitle("c#t (#mum)");
+       hmultiplycT->GetYaxis()->SetTitleOffset(1.5);
+       cmultiplycT->SetLeftMargin(0.15);
+       cmultiplycT->SetRightMargin(0.05);
+       cmultiplycT->cd();
+       hmultiplycT->SetFillStyle(3004);
+       hmultiplycT->SetFillColor(38);
+       hmultiplycT->Draw("hist err");
+       hMCcT1D->Draw("hist same err");
+
+       /*
+       TFile* file_histo = new TFile("fileHisto_180004.root","RECREATE");
+       hMCpt1D->Write("hMCpt1D");
+       hRECpt1D->Write("hRECpt1D");
+       hMCy1D->Write("hMCy1D");
+       hRECy1D->Write("hRECy1D");
+       hMCcTS1D->Write("hMCcTS1D");
+       hRECcTS1D->Write("hRECcTS1D");
+       hMCptPi1D->Write("hMCptPi1D");
+       hRECptPi1D->Write("hRECptPi1D");
+       hMCptK1D->Write("hMCptK1D");
+       hRECptK1D->Write("hRECptK1D");
+       hMCcT1D->Write("hMCcT1D");
+       hRECcT1D->Write("hRECcT1D");
+       hpteffCF->Write("hpteffCF");
+       hyeffCF->Write("hyeffCF");
+       hcTSeffCF->Write("hcTSeffCF");
+       hptPieffCF->Write("hptPieffCF");
+       hptKeffCF->Write("hptKeffCF");
+       hcTeffCF->Write("hcTeffCF");
+       hmultiplypt->Write("");
+       hmultiplyy->Write("");
+       hmultiplycTS->Write("");
+       hmultiplyptPi->Write("");
+       hmultiplyptK->Write("");
+       hmultiplycT->Write("");
+       file_histo->Close();
+       */
+
+       Int_t nsliceVars = 5;
+       Int_t sliceVars[5];
+       sliceVars[0] = 1; 
+       sliceVars[1] = 2; 
+       sliceVars[2] = 3; 
+       sliceVars[3] = 4; 
+       sliceVars[4] = 5; 
+       Double_t varMin[6] = {0.5, -2, -1, 0., 0., 0. };
+       Double_t varMax[6] = {1.5, 2, 1, 10., 10., 500.};
+       AliCFContainer* data_sliced = data->MakeSlice(nsliceVars, sliceVars, varMin, varMax);
+       cout << " the container now has " << data_sliced->GetNVar() << " dimensions " << endl;
+       
+}