New task for D* efficiencies with CORRFW (Alessandro)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 18 Feb 2010 14:11:47 +0000 (14:11 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 18 Feb 2010 14:11:47 +0000 (14:11 +0000)
PWG3/PWG3vertexingHFLinkDef.h
PWG3/libPWG3vertexingHF.pkg
PWG3/vertexingHF/AddTaskCFDStar.C [new file with mode: 0644]
PWG3/vertexingHF/AliCFTaskForDStarAnalysis.cxx [new file with mode: 0644]
PWG3/vertexingHF/AliCFTaskForDStarAnalysis.h [new file with mode: 0644]

index 239477b..3224fdc 100644 (file)
@@ -20,6 +20,7 @@
 #pragma link C++ class AliAnalysisTaskSED0Mass+;
 #pragma link C++ class AliAnalysisTaskSECharmFraction+;
 #pragma link C++ class AliCFHeavyFlavourTaskMultiVarMultiStep+;
+#pragma link C++ class AliCFTaskForDStarAnalysis+;
 #pragma link C++ class AliAnalysisTaskSEDStar+;
 #pragma link C++ class AliAnalysisTaskSEDStarJets+;
 #pragma link C++ class AliMultiDimVector+;
index 2e6e33b..b746b41 100644 (file)
@@ -14,6 +14,7 @@ SRCS:=   vertexingHF/AliAODRecoDecayHF.cxx \
        vertexingHF/AliAnalysisTaskSED0Mass.cxx \
        vertexingHF/AliAnalysisTaskSECharmFraction.cxx \
        vertexingHF/AliCFHeavyFlavourTaskMultiVarMultiStep.cxx \
+       vertexingHF/AliCFTaskForDStarAnalysis.cxx \
        vertexingHF/AliAnalysisTaskSEDStar.cxx \
        vertexingHF/AliAnalysisTaskSEDStarJets.cxx \
        vertexingHF/AliMultiDimVector.cxx vertexingHF/AliSignificanceCalculator.cxx \
diff --git a/PWG3/vertexingHF/AddTaskCFDStar.C b/PWG3/vertexingHF/AddTaskCFDStar.C
new file mode 100644 (file)
index 0000000..4be370e
--- /dev/null
@@ -0,0 +1,382 @@
+//DEFINITION OF A FEW CONSTANTS
+//
+// binning method from C.Zampolli
+//
+// general
+const Double_t ymin  = -2.1 ;
+const Double_t ymax  =  2.1 ;
+//soft pion
+const Double_t ptmin_0_1 =  0.0 ;
+const Double_t ptmax_0_1 =  1.0 ;
+const Double_t ptmin_1_2 =  1.0 ;
+const Double_t ptmax_1_2 =  2.0 ;
+const Double_t ptmin_2_10 =  2.0 ;
+const Double_t ptmax_2_10 =  15.0 ;
+//D0 and D0 prongs
+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 =  20.0 ;
+const Double_t cosmin = -1.05;
+const Double_t cosmax =  1.05;
+const Double_t cTmin = 0;  // micron
+const Double_t cTmax = 500;  // micron
+const Double_t dcamin = 0;  // micron
+const Double_t dcamax = 500;  // micron
+const Double_t d0min = -1000;  // micron
+const Double_t d0max = 1000;  // micron
+const Double_t d0xd0min = -100000;  // micron
+const Double_t d0xd0max = 100000;  // micron
+const Double_t phimin = 0.0;    
+const Int_t    mintrackrefsTPC = 2 ;
+const Int_t    mintrackrefsITS = 3 ;
+const Int_t    charge  = 1 ; 
+const Int_t    minclustersTPC = 50 ;
+// cuts
+const Double_t ptmin = 0.05;
+const Double_t ptmax = 9999.;
+const Double_t etamin = -0.9;
+const Double_t etamax = 0.9;
+const Double_t zmin = -15;
+const Double_t zmax = 15;
+const Int_t    minITSClusters = 3;
+const Int_t    minITSClustersSoft = 2;
+//----------------------------------------------------
+
+AliCFTaskForDStarAnalysis *AddTaskCFDStar()
+{
+
+       //CONTAINER DEFINITION
+       Info("AliCFTaskForDStarAnalysis","SETUP CONTAINER");
+       //the sensitive variables, their indices
+       UInt_t ipt = 0;
+       UInt_t iy  = 1;
+       UInt_t icosThetaStar  = 2;
+       UInt_t ipTpi  = 3;
+       UInt_t ipTD0  = 4;
+       UInt_t icT    = 5;
+       UInt_t idca   = 6;
+       UInt_t id0pi  = 7;
+       UInt_t id0K   = 8;
+       UInt_t id0xd0     = 9;
+       UInt_t ipointing  = 10;
+       UInt_t iphi  = 11;
+       UInt_t iz    = 12;
+        UInt_t ipTD0pi = 13;
+        UInt_t ipTD0K  = 14;
+
+       const Double_t phimax = 2*TMath::Pi();
+
+       //Setting up the container grid... 
+       UInt_t nstep = 8; //number of selection steps
+       const Int_t nvar   = 15 ; //number of variables on the grid:pt, y, cosThetaStar, pTpi, pTk, cT, dca, d0pi, d0K, d0xd0, cosPointingAngle, phi 
+       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  = 2 ; //bins in pt from 8 to 10 GeV
+       const Int_t nbin1  = 30 ; //bins in y
+       const Int_t nbin2  = 30 ; //bins in cosThetaStar 
+        // soft pion and D0 from D*
+       const Int_t nbin3_0_1  = 8 ; //bins in ptPi from 0 to 4 GeV
+       const Int_t nbin3_1_2  = 1 ; //bins in ptPi from 4 to 8 GeV
+       const Int_t nbin3_2_10 = 1 ; //bins in ptPi from 8 to 10 GeV
+       const Int_t nbin4_0_4  = 8 ; //bins in ptD0 from 0 to 4 GeV
+       const Int_t nbin4_4_8  = 3 ; //bins in ptD0 from 4 to 8 GeV
+       const Int_t nbin4_8_10 = 1 ; //bins in ptD0 from 8 to 10 GeV
+        // D0 prongs - cutting variables
+       const Int_t nbin5  = 20 ; //bins in cT
+       const Int_t nbin6  = 20 ; //bins in dca
+       const Int_t nbin7  = 100 ; //bins in d0pi
+       const Int_t nbin8  = 100 ; //bins in d0K
+       const Int_t nbin9  = 80 ; //bins in d0xd0
+       const Int_t nbin10  = 100 ; //bins in cosPointingAngle
+       const Int_t nbin11  = 15 ; //bins in Phi
+               const Int_t nbin12  = 60 ; //bins in z vertex
+        // D0 prongs pt and phi
+       const Int_t nbin5_0_4   = 8 ; //bins in ptPi from 0 to 4 GeV
+       const Int_t nbin5_4_8   = 4 ; //bins in ptPi from 4 to 8 GeV
+       const Int_t nbin5_8_10  = 8 ; //bins in ptPi from 8 to 10 GeV
+       const Int_t nbin6_0_4   = 8 ; //bins in ptk from 0 to 4 GeV
+       const Int_t nbin6_4_8   = 4 ; //bins in ptk from 4 to 8 GeV
+       const Int_t nbin6_8_10  = 8 ; //bins in ptk from 8 to 10 GeV
+        
+       //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_1+nbin3_1_2+nbin3_2_10;
+       iBin[4]=nbin4_0_4+nbin4_4_8+nbin4_8_10;
+       iBin[5]=nbin5;
+       iBin[6]=nbin6;
+       iBin[7]=nbin7;
+       iBin[8]=nbin8;
+       iBin[9]=nbin9;
+       iBin[10]=nbin10;
+       iBin[11]=nbin11;
+       iBin[12]=nbin12;
+        iBin[13]=nbin5_0_4+nbin5_4_8+nbin5_8_10;
+       iBin[14]=nbin6_0_4+nbin6_4_8+nbin6_8_10;
+       
+       //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];
+       Double_t *binLim6 = new Double_t[iBin[6]+1];
+       Double_t *binLim7 = new Double_t[iBin[7]+1];
+       Double_t *binLim8 = new Double_t[iBin[8]+1];
+       Double_t *binLim9 = new Double_t[iBin[9]+1];
+       Double_t *binLim10 = new Double_t[iBin[10]+1];
+       Double_t *binLim11 = new Double_t[iBin[11]+1];
+       Double_t *binLim12 = new Double_t[iBin[12]+1];
+       Double_t *binLim13 = new Double_t[iBin[13]+1];
+       Double_t *binLim14 = new Double_t[iBin[14]+1];
+
+       // checking limits
+       if (ptmax_0_4 != ptmin_4_8) {
+               Error("AliCFTaskForDStarAnalysis","max lim 1st range != min lim 2nd range, please check!");
+       }
+       if (ptmax_4_8 != ptmin_8_10) {
+               Error("AliCFTaskForDStarAnalysis","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("AliCFDStar","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("AliCFDStar","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 ;
+
+       // Soft ptPi ---------------------------------------------------------------------------------
+       for(Int_t i=0; i<=nbin3_0_1; i++) binLim3[i]=(Double_t)ptmin_0_1 + (ptmax_0_1-ptmin_0_1)/nbin3_0_1*(Double_t)i ; 
+       if (binLim3[nbin3_0_1] != ptmin_1_2)  {
+               Error("AliCFDStar","Calculated bin lim for ptPi - 1st range - differs from expected!");
+       }
+       for(Int_t i=0; i<=nbin3_1_2; i++) binLim3[i+nbin3_0_1]=(Double_t)ptmin_1_2 + (ptmax_1_2-ptmin_1_2)/nbin3_1_2*(Double_t)i ; 
+       if (binLim3[nbin3_0_1+nbin3_1_2] != ptmin_2_10)  {
+               Error("AliCFDStar","Calculated bin lim for ptPi - 2nd range - differs from expected!\n");
+       }
+       for(Int_t i=0; i<=nbin3_2_10; i++) binLim3[i+nbin3_0_1+nbin3_1_2]=(Double_t)ptmin_2_10 + (ptmax_2_10-ptmin_2_10)/nbin3_2_10*(Double_t)i ; 
+
+       // ptD0 --------------------------------------------------------------------------------------
+       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("AliCFDStar","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("AliCFDStar","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 ; 
+        // D0 ptPi --------------------------------------------------------------------------------------------------------
+        for(Int_t i=0; i<=nbin5_0_4; i++) binLim13[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin5_0_4*(Double_t)i ; 
+        if (binLim13[nbin5_0_4] != ptmin_4_8)  {
+                Error("AliCFDStar","Calculated bin lim for ptPi - 1st range - differs from expected!");
+        }
+        for(Int_t i=0; i<=nbin5_4_8; i++) binLim13[i+nbin5_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin5_4_8*(Double_t)i ; 
+        if (binLim13[nbin5_0_4+nbin5_4_8] != ptmin_8_10)  {
+                Error("AliCFDStar","Calculated bin lim for ptPi - 2nd range - differs from expected!\n");
+        }
+        for(Int_t i=0; i<=nbin5_8_10; i++) binLim13[i+nbin5_0_4+nbin5_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin5_8_10*(Double_t)i ; 
+
+        // D0 ptK ----------------------------------------------------------------------------------------------------------
+        for(Int_t i=0; i<=nbin6_0_4; i++) binLim14[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin6_0_4*(Double_t)i ; 
+        if (binLim14[nbin6_0_4] != ptmin_4_8)  {
+                Error("AliCFDStar","Calculated bin lim for ptKa - 1st range - differs from expected!");
+        }
+        for(Int_t i=0; i<=nbin6_4_8; i++) binLim14[i+nbin6_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin6_4_8*(Double_t)i ; 
+        if (binLim14[nbin6_0_4+nbin6_4_8] != ptmin_8_10)  {
+                Error("AliCFDStar","Calculated bin lim for ptKa - 2nd range - differs from expected!\n");
+        }
+        for(Int_t i=0; i<=nbin6_8_10; i++) binLim14[i+nbin6_0_4+nbin6_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin6_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 ;
+
+       // dca
+       for(Int_t i=0; i<=nbin6; i++) binLim6[i]=(Double_t)dcamin  + (dcamax-dcamin)  /nbin6*(Double_t)i ;
+
+       // d0pi
+       for(Int_t i=0; i<=nbin7; i++) binLim7[i]=(Double_t)d0min  + (d0max-d0min)  /nbin7*(Double_t)i ;
+
+       // d0K
+       for(Int_t i=0; i<=nbin8; i++) binLim8[i]=(Double_t)d0min  + (d0max-d0min)  /nbin8*(Double_t)i ;
+
+       // d0xd0
+       for(Int_t i=0; i<=nbin9; i++) binLim9[i]=(Double_t)d0xd0min  + (d0xd0max-d0xd0min)  /nbin9*(Double_t)i ;
+
+       // cosPointingAngle
+       for(Int_t i=0; i<=nbin10; i++) binLim10[i]=(Double_t)cosmin  + (cosmax-cosmin)  /nbin10*(Double_t)i ;
+
+       // Phi
+       for(Int_t i=0; i<=nbin11; i++) binLim11[i]=(Double_t)phimin  + (phimax-phimin)  /nbin11*(Double_t)i ;
+
+       // z Primary Vertex
+       for(Int_t i=0; i<=nbin12; i++) {
+               binLim12[i]=(Double_t)zmin  + (zmax-zmin)  /nbin12*(Double_t)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(ipTD0,binLim4);
+       container -> SetBinLimits(icT,binLim5);
+       container -> SetBinLimits(idca,binLim6);
+       container -> SetBinLimits(id0pi,binLim7);
+       container -> SetBinLimits(id0K,binLim8);
+       container -> SetBinLimits(id0xd0,binLim9);
+       container -> SetBinLimits(ipointing,binLim10);
+       container -> SetBinLimits(iphi,binLim11);
+       container -> SetBinLimits(iz,binLim12);
+       container -> SetBinLimits(ipTD0pi,binLim13);
+       container -> SetBinLimits(ipTD0K,binLim14);
+       
+       //CREATE THE  CUTS -----------------------------------------------
+       
+       // Gen-Level kinematic cuts
+       AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
+       
+       //Particle-Level cuts:  
+       AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
+       mcGenCuts->SetRequirePdgCode(413, kTRUE);  // kTRUE set in order to include D*_bar
+       mcGenCuts->SetAODMC(1); //special flag for reading MC in AOD tree (important)
+       
+       // Acceptance cuts:
+       AliCFAcceptanceCuts* accCuts = new AliCFAcceptanceCuts("accCuts", "Acceptance cuts");
+       AliCFTrackKineCuts *kineAccCuts = new AliCFTrackKineCuts("kineAccCuts","Kine-Acceptance cuts");
+       kineAccCuts->SetPtRange(ptmin,ptmax);
+       kineAccCuts->SetEtaRange(etamin,etamax);
+
+       // Rec-Level kinematic cuts
+       AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");
+       
+       AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
+       
+       AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
+       
+       printf("CREATE MC KINE CUTS\n");
+       TObjArray* mcList = new TObjArray(0) ;
+       mcList->AddLast(mcKineCuts);
+       mcList->AddLast(mcGenCuts);
+       
+       printf("CREATE ACCEPTANCE CUTS\n");
+       TObjArray* accList = new TObjArray(0) ;
+       accList->AddLast(kineAccCuts);
+
+       printf("CREATE RECONSTRUCTION CUTS\n");
+       TObjArray* recList = new TObjArray(0) ;   // not used!! 
+       recList->AddLast(recKineCuts);
+       recList->AddLast(recQualityCuts);
+       recList->AddLast(recIsPrimaryCuts);
+       
+       TObjArray* emptyList = new TObjArray(0);
+
+       //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 , accList); // Acceptance 
+       man->SetParticleCutsList(2 , emptyList); // Vertex 
+       man->SetParticleCutsList(3 , emptyList); // Refit 
+       man->SetParticleCutsList(4 , emptyList); // AOD
+       man->SetParticleCutsList(5 , emptyList); // AOD in Acceptance
+       man->SetParticleCutsList(6 , emptyList); // AOD with required n. of ITS clusters
+       man->SetParticleCutsList(7 , emptyList); // AOD Reco cuts
+       
+       // Get the pointer to the existing analysis manager via the static access method.
+       //==============================================================================
+       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+       if (!mgr) {
+         ::Error("AddTaskCompareHF", "No analysis manager to connect to.");
+         return NULL;
+       }   
+       //CREATE THE TASK
+       printf("CREATE TASK\n");
+       // create the task
+       AliCFTaskForDStarAnalysis *task = new AliCFTaskForDStarAnalysis("AliCFTaskForDStarAnalysis");
+       task->SetMinITSClusters(minITSClusters);
+       task->SetMinITSClustersSoft(minITSClustersSoft);
+       task->SetCFManager(man); //here is set the CF manager
+       
+        Bool_t AcceptanceUnf = kTRUE; // unfold at acceptance level, otherwise D* cuts
+        Int_t thnDim[4];
+        
+        //first half  : reconstructed 
+        //second half : MC
+        thnDim[0] = iBin[0];
+        thnDim[2] = iBin[0];
+        thnDim[1] = iBin[1];
+        thnDim[3] = iBin[1];
+
+        THnSparseD* correlation = new THnSparseD("correlation","THnSparse with correlations",4,thnDim);
+        Double_t** binEdges = new Double_t[2];
+
+        // set bin limits
+
+        binEdges[0]= binLim0;
+        binEdges[1]= binLim1;
+
+        correlation->SetBinEdges(0,binEdges[0]);
+        correlation->SetBinEdges(2,binEdges[0]);
+
+        correlation->SetBinEdges(1,binEdges[1]);
+        correlation->SetBinEdges(3,binEdges[1]);
+
+        correlation->Sumw2();
+  
+        // correlation matrix ready
+        //------------------------------------------------//
+
+        task->SetCorrelationMatrix(correlation); // correlation matrix for unfolding
+       
+       // Create and connect containers for input/output
+       
+       // ------ input data ------
+       AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();
+       
+       // ----- output data -----
+       
+       TString outputfile = AliAnalysisManager::GetCommonFileName();
+       outputfile += ":PWG3_D2H_CFtaskDStar";
+
+       //now comes user's output objects :
+       // output TH1I for event counting
+       AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("CFDSchist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+       // output Correction Framework Container (for acceptance & efficiency calculations)
+       AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("CFDSccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+       // Unfolding - correlation matrix
+        AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("CFDScorr0", THnSparseD::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+
+       mgr->AddTask(task);
+       
+       mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+       mgr->ConnectOutput(task,1,coutput1);
+       mgr->ConnectOutput(task,2,coutput2);
+        mgr->ConnectOutput(task,3,coutput3);
+
+       return task;
+}
+
diff --git a/PWG3/vertexingHF/AliCFTaskForDStarAnalysis.cxx b/PWG3/vertexingHF/AliCFTaskForDStarAnalysis.cxx
new file mode 100644 (file)
index 0000000..2874fa2
--- /dev/null
@@ -0,0 +1,753 @@
+/**************************************************************************
+ * 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 DStar corrections: 
+//   
+//   The D0 cutting varibles position in the container and container 
+//   binning method from a C.Zampolli example
+//   In this way a simple comparison between D0 and D0 from D* is possible. 
+//
+//-----------------------------------------------------------------------
+// Author : A.Grelli, Utrecht University
+//
+//         a.grelli@uu.nl
+//-----------------------------------------------------------------------
+#include <TCanvas.h>
+#include <TParticle.h>
+#include <TDatabasePDG.h>
+#include <TH1I.h>
+#include <TStyle.h>
+#include <TFile.h>
+
+#include "AliCFTaskForDStarAnalysis.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliCFManager.h"
+#include "AliCFContainer.h"
+#include "AliLog.h"
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODRecoDecay.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "AliESDtrack.h"
+#include "TChain.h"
+#include "THnSparse.h"
+#include "TH2D.h"
+
+//__________________________________________________________________________
+AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis() :
+  AliAnalysisTaskSE(),
+  fCFManager(0x0),
+  fHistEventsProcessed(0x0),
+  fCorrelation(0x0),
+  fCountRecoDStarSel(0),
+  fEvents(0),
+  fMinITSClusters(5),
+  fMinITSClustersSoft(4),
+  fAcceptanceUnf(kTRUE)
+{
+  //
+  //Default ctor
+  //
+}
+//___________________________________________________________________________
+AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis(const Char_t* name) :
+  AliAnalysisTaskSE(name),
+  fCFManager(0x0),
+  fHistEventsProcessed(0x0),
+  fCorrelation(0x0),
+  fCountRecoDStarSel(0),
+  fEvents(0),
+  fMinITSClusters(5),
+  fMinITSClustersSoft(4),
+  fAcceptanceUnf(kTRUE)
+{
+  //
+  // Constructor. Initialization of Inputs and Outputs
+  //
+  Info("AliCFTaskForDStarAnalysis","Calling Constructor");
+  
+  DefineOutput(1,TH1I::Class());
+  DefineOutput(2,AliCFContainer::Class());
+  DefineOutput(3,THnSparseD::Class());
+}
+
+//___________________________________________________________________________
+AliCFTaskForDStarAnalysis& AliCFTaskForDStarAnalysis::operator=(const AliCFTaskForDStarAnalysis& c) 
+{
+  //
+  // Assignment operator
+  //
+  if (this!=&c) {
+    AliAnalysisTaskSE::operator=(c) ;
+    fCFManager  = c.fCFManager;
+    fHistEventsProcessed = c.fHistEventsProcessed;
+  }
+  return *this;
+}
+
+//___________________________________________________________________________
+AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis(const AliCFTaskForDStarAnalysis& c) :
+  AliAnalysisTaskSE(c),
+  fCFManager(c.fCFManager),
+  fHistEventsProcessed(c.fHistEventsProcessed),
+  fCorrelation(c.fCorrelation),
+  fCountRecoDStarSel(c.fCountRecoDStarSel),
+  fEvents(c.fEvents),
+  fMinITSClusters(c.fMinITSClusters),
+  fMinITSClustersSoft(c.fMinITSClustersSoft),
+  fAcceptanceUnf(c.fAcceptanceUnf)
+{
+  //
+  // Copy Constructor
+  //
+}
+
+//___________________________________________________________________________
+AliCFTaskForDStarAnalysis::~AliCFTaskForDStarAnalysis() {
+  //
+  //destructor
+  //
+  if (fCFManager)           delete fCFManager ;
+  if (fHistEventsProcessed) delete fHistEventsProcessed ;
+  if (fCorrelation)         delete fCorrelation ;
+}
+
+//_________________________________________________
+void AliCFTaskForDStarAnalysis::UserExec(Option_t *)
+{
+  //
+  // Main loop function
+  //
+  
+  if (!fInputEvent) {
+    Error("UserExec","NO EVENT FOUND!");
+    return;
+  }
+  
+  fEvents++;
+
+  if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
+  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
+  
+  TClonesArray *arrayDStartoD0pi=0;
+  
+  if(!aodEvent && AODEvent() && IsStandardAOD()) {
+    // In case there is an AOD handler writing a standard AOD, use the AOD 
+    // event in memory rather than the input (ESD) event.    
+    aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
+    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
+    // have to taken from the AOD event hold by the AliAODExtension
+    AliAODHandler* aodHandler = (AliAODHandler*) 
+      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+    if(aodHandler->GetExtensions()) {
+      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
+      AliAODEvent *aodFromExt = ext->GetAOD();
+      arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
+    }
+  } else {
+    arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
+  }
+  
+  if (!arrayDStartoD0pi) {
+    AliError("Could not find array of HF vertices");
+    return;
+  }
+  
+  fCFManager->SetRecEventInfo(aodEvent);
+  fCFManager->SetMCEventInfo(aodEvent);
+  
+  // event selection
+  Double_t containerInput[14] ;
+  Double_t containerInputMC[14] ;
+  
+  //loop on the MC event
+  
+  TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+  if (!mcArray) {
+    AliError("Could not find Monte-Carlo in AOD");
+    return;
+  }
+  
+  AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+  if (!mcHeader) {
+    AliError("Could not find MC Header in AOD");
+    return;
+  }
+  
+  // AOD primary vertex
+  AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
+  Double_t zPrimVertex = vtx1->GetZ();
+  Double_t zMCVertex = mcHeader->GetVtxZ();
+  Bool_t vtxFlag = kTRUE;
+  TString title=vtx1->GetTitle();
+  if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
+  
+  for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { 
+    AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
+    if (!mcPart) {
+      AliWarning("MC 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[9] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.};
+    
+    if (GetDStarMCParticle(mcPart, mcArray, vectorMC)){
+      
+      containerInputMC[0] = vectorMC[0];
+      containerInputMC[1] = vectorMC[1] ;
+      containerInputMC[2] = vectorMC[2] ;
+      containerInputMC[3] = vectorMC[3] ;
+      containerInputMC[4] = vectorMC[4] ;
+      containerInputMC[5] = vectorMC[5] ; 
+      containerInputMC[6] = 0.;      
+      containerInputMC[7] = 0.;          
+      containerInputMC[8] = 0.;
+      containerInputMC[9] = -100000.; 
+      containerInputMC[10] = 1.01; 
+      containerInputMC[11] = vectorMC[6];  
+      containerInputMC[12] = zMCVertex;     // z vertex
+      containerInputMC[13] = vectorMC[7];   // pt D0 pion
+      containerInputMC[14] = vectorMC[8];   // pt D0 kaon  
+      
+      fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated);
+      
+      // check the MC-Acceptance level cuts
+      // since standard CF functions are not applicable, using Kine Cuts on daughters
+      
+      Int_t daughter0 = mcPart->GetDaughter(0);
+      Int_t daughter1 = mcPart->GetDaughter(1);
+      
+      AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
+      
+      //D0
+      AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
+      // Soft Pion
+      AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
+      
+      // Acceptance variables for the soft pion
+      Double_t eta1 = mcPartDaughter1->Eta();
+      Double_t pt1  = mcPartDaughter1->Pt(); 
+      
+      Int_t daughD00 = 0;
+      Int_t daughD01 = 0;
+      
+      // Just to be sure to take the right particles    
+      if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){
+       daughD00 = mcPartDaughter0->GetDaughter(0);
+       daughD01 = mcPartDaughter0->GetDaughter(1);
+      }else{
+       daughD00 = mcPartDaughter1->GetDaughter(0);
+       daughD01 = mcPartDaughter1->GetDaughter(1);
+      }
+      
+      // the D0 daughters
+      AliAODMCParticle* mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00));
+      AliAODMCParticle* mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01));
+      
+      if (!mcD0PartDaughter0 || !mcD0PartDaughter1) {
+       AliWarning("At least one D0 Daughter Particle not found in tree, but it should be, this check was already done..."); 
+      }
+      
+      // D0 daughters - needed for acceptance
+      Double_t theD0pt0 =  mcD0PartDaughter0->Pt();
+      Double_t theD0pt1 =  mcD0PartDaughter1->Pt();
+      Double_t theD0eta0 = mcD0PartDaughter0->Eta();
+      Double_t theD0eta1 = mcD0PartDaughter1->Eta();
+      
+      // ACCEPTANCE REQUESTS ---------
+      
+      // soft pion 
+      Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 >= 0.05);
+      // Do daughters
+      Bool_t D0daught0inAcceptance = (TMath::Abs(theD0eta0) <= 0.9 && theD0pt0 >= 0.1); 
+      Bool_t D0daught1inAcceptance = (TMath::Abs(theD0eta1) <= 0.9 && theD0pt1 >= 0.1);
+      
+      if (daught1inAcceptance && D0daught0inAcceptance && D0daught1inAcceptance) {
+       
+       AliDebug(2, "D* Daughter particles in acceptance");
+       
+       fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
+       
+       // check on the vertex
+       if (vtxFlag){
+         // filling the container if the vertex is ok
+         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
+         
+         Bool_t refitFlag = kTRUE;
+         for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
+           AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
+           
+            // refit only for D0 daughters
+           if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) {
+             if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
+               refitFlag = kFALSE;
+             }
+           }
+         } 
+         if (refitFlag){ // refit
+      
+           fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
+
+         } // end of refit       
+       } // end of vertex
+      }        //end of acceptance             
+    } // end of MC D*
+    else {
+      AliDebug(3,"Problems in filling the container");
+      continue;
+    }
+  } // end of MC loop
+  
+  //rec
+  AliDebug(2, Form("Found %d D* candidates",arrayDStartoD0pi->GetEntriesFast()));
+  
+  //D* and D0 prongs needed to MatchToMC method
+  Int_t pdgDgDStartoD0pi[2]={421,211};
+  Int_t pdgDgD0toKpi[2]={321,211};
+  
+  for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
+    
+    // D* candidates
+    AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
+    
+    // D0 from the reco cascade
+    AliAODRecoDecayHF2Prong* D0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
+    Bool_t unsetvtx=kFALSE;
+    
+    // needed for pointing angle
+    if(!D0particle->GetOwnPrimaryVtx()) {
+      D0particle->SetOwnPrimaryVtx(vtx1);
+      unsetvtx=kTRUE;
+    }
+    
+    // find associated MC particle for D* ->D0toKpi
+    Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray); 
+    
+    // find D0->Kpi ... needed in the following
+    Int_t mcD0Label = D0particle->MatchToMC(421,mcArray,2,pdgDgD0toKpi); 
+    
+    if (mcLabel == -1 || mcD0Label ==-1) 
+      {
+       AliDebug(2,"No MC particle found");
+      }
+    else {
+      
+      // the D* and the D0 in MC
+      AliAODMCParticle* mcVtxHFDStar = (AliAODMCParticle*)mcArray->At(mcLabel);
+      AliAODMCParticle* mcVtxHFD0Kpi = (AliAODMCParticle*)mcArray->At(mcD0Label);
+      
+      if (!mcVtxHFD0Kpi || !mcVtxHFDStar) {
+       AliWarning("Could not find associated MC D0 and/or D* in AOD MC tree");
+       continue;
+      }
+      
+      // soft pion
+      AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor(); 
+
+      //D0tokpi
+      AliAODTrack *track0 = (AliAODTrack*)D0particle->GetDaughter(0);
+      AliAODTrack *track1 = (AliAODTrack*)D0particle->GetDaughter(1);
+      
+      // check if associated MC v0 passes the cuts
+      if (!fCFManager->CheckParticleCuts(0 ,mcVtxHFDStar)) { 
+       AliDebug(2, "Skipping the particles due to cuts");
+       continue; 
+      }
+      
+      // fill the container...
+      Double_t pt = TMath::Sqrt(TMath::Power(D0particle->Pt(),2)+TMath::Power(track2->Pt(),2));
+      Double_t rapidity = dstarD0pi->YDstar();
+      Double_t cosThetaStar = 9999.;
+      Double_t pTpi  = 0.;
+      Double_t pTK   = 0.;
+      Double_t dca   = (D0particle->GetDCA())*1E4;
+      Double_t d0pi  = 0.;
+      Double_t d0K   = 0.;
+      Double_t d0xd0 = (D0particle->Prodd0d0())*1E8;
+      Double_t cosPointingAngle = D0particle->CosPointingAngle();
+      Double_t phi = D0particle->Phi();
+           
+      // Select D0 cutting variables
+      Int_t pdgCode = mcVtxHFD0Kpi->GetPdgCode();
+
+      // D0 related variables
+      if (pdgCode > 0){
+
+       cosThetaStar = D0particle->CosThetaStarD0();
+       pTpi    = D0particle->PtProng(0);
+       pTK     = D0particle->PtProng(1);
+       d0pi    = (D0particle->Getd0Prong(0))*1E4;
+       d0K     = (D0particle->Getd0Prong(1))*1E4;
+   
+      }
+      else {
+
+       cosThetaStar = D0particle->CosThetaStarD0bar();
+       pTpi    = D0particle->PtProng(1);
+       pTK     = D0particle->PtProng(0);
+       d0pi    = (D0particle->Getd0Prong(1))*1E4;
+       d0K     = (D0particle->Getd0Prong(0))*1E4;
+   
+      }
+      
+      // ct of the D0 from D*
+      Double_t cT = D0particle->CtD0();
+      
+      containerInput[0]  = pt;
+      containerInput[1]  = rapidity;
+      containerInput[2]  = cosThetaStar;
+      containerInput[3]  = track2->Pt();
+      containerInput[4]  = D0particle->Pt();
+      containerInput[5]  = cT*1.E4;  // in micron
+      containerInput[6]  = dca;  // in micron
+      containerInput[7]  = d0pi;  // in micron
+      containerInput[8]  = d0K;  // in micron
+      containerInput[9]  = d0xd0;  // in micron^2
+      containerInput[10] = cosPointingAngle;  // in micron
+      containerInput[11] = phi;  
+      containerInput[12] = zPrimVertex;    // z of reconstructed of primary vertex
+      containerInput[13] = pTpi;  // D0 pion
+      containerInput[14] = pTK;   // D0 kaon  
+
+      fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
+      
+      // refit in ITS and TPC for D0 daughters
+      if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
+       continue;
+      }
+
+      // reft in ITS for soft pion
+      if((!(track2->GetStatus()&AliESDtrack::kITSrefit))) {
+       continue;
+      }
+
+      // cut in acceptance for the soft pion and for the D0 daughters      
+      Bool_t acceptanceProng0 = (TMath::Abs(D0particle->EtaProng(0))<= 0.9 && D0particle->PtProng(0) >= 0.1);
+      Bool_t acceptanceProng1 = (TMath::Abs(D0particle->EtaProng(1))<= 0.9 && D0particle->PtProng(1) >= 0.1);
+      
+      // soft pion acceptance ... is it fine 0.9?????
+      Bool_t acceptanceProng2 = (TMath::Abs(track2->Eta())<= 0.9 && track2->Pt() >= 0.05);
+      
+      if (acceptanceProng0 && acceptanceProng1 && acceptanceProng2) {
+       AliDebug(2,"D* reco daughters in acceptance");
+       fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
+       
+       if(fAcceptanceUnf){       
+         Double_t fill[4]; //fill response matrix
+         
+         // dimensions 0&1 : pt,eta (Rec)              
+         fill[0] = pt ;
+         fill[1] = rapidity;           
+         // dimensions 2&3 : pt,eta (MC)                                       
+         fill[2] =  mcVtxHFDStar->Pt();
+         fill[3] =  mcVtxHFDStar->Y();   
+         fCorrelation->Fill(fill);             
+       }  
+       
+       // cut on the min n. of clusters in ITS for the D0 and soft pion
+       Int_t ncls0=0,ncls1=0,ncls2=0;
+       for(Int_t l=0;l<6;l++) {
+         if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
+         if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
+         if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
+       }       
+       // see AddTask for soft pion ITS clusters request
+       AliDebug(2, Form("n clusters = %d", ncls0));
+
+       if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>= fMinITSClustersSoft) {
+         fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
+         
+         // D0 cuts optimized for D* analysis
+         Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027}; 
+
+          // needed for cuts     
+          Double_t theD0pt = D0particle->Pt();
+
+         if (theD0pt <= 1){ // first bin not optimized
+           cuts[0] = 400;
+           cuts[1] = 0.8;
+           cuts[2] = 0.21;
+           cuts[3] = 500;
+           cuts[4] = 500;
+           cuts[5] = -20000;
+           cuts[6] = 0.6;  
+         }
+         else if (theD0pt > 1 && theD0pt <= 2){
+           cuts[0] = 200; 
+           cuts[1] = 0.7; 
+           cuts[2] = 0.8; 
+           cuts[3] = 210;
+           cuts[4] = 210;
+           cuts[5] = -20000;
+           cuts[6] = 0.9;  
+         }
+         else if (theD0pt > 2 && theD0pt <= 3){
+           cuts[0] = 400;
+           cuts[1] = 0.8; 
+           cuts[2] = 0.8;
+           cuts[3] = 420;
+           cuts[4] = 350; 
+           cuts[5] = -8500;
+           cuts[6] = 0.9;   
+         }
+         else if (theD0pt > 3 && theD0pt <= 5){
+           cuts[0] = 160;  
+           cuts[1] = 1.0; 
+           cuts[2] = 1.2;  
+           cuts[3] = 560; 
+           cuts[4] = 420; 
+           cuts[5] = -8500;
+           cuts[6] = 0.9;  
+         }
+         else if (theD0pt > 5){
+           cuts[0] = 800;
+           cuts[1] = 1.0;
+           cuts[2] = 1.2; 
+           cuts[3] = 700;
+           cuts[4] = 700; 
+           cuts[5] = 10000;
+           cuts[6] = 0.9;  
+         }
+         if (dca < cuts[0] 
+             && TMath::Abs(cosThetaStar) < cuts[1]  
+             && pTpi > cuts[2] 
+             && pTK > cuts[2]  
+             && TMath::Abs(d0pi) < cuts[3] 
+             && TMath::Abs(d0K) < cuts[4]  
+             && d0xd0 < cuts[5] 
+             && cosPointingAngle > cuts[6]
+           ){
+           
+           AliDebug(2,"Particle passed D* selection cuts");
+           fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoCuts) ;   
+           
+           if(!fAcceptanceUnf){ // unfolding
+             
+             Double_t fill[4]; //fill response matrix
+             
+             // dimensions 0&1 : pt,eta (Rec)              
+             fill[0] = pt ;
+             fill[1] = rapidity;                   
+             // dimensions 2&3 : pt,eta (MC)               
+             fill[2] =  mcVtxHFDStar->Pt();
+             fill[3] =  mcVtxHFDStar->Y();
+             
+             fCorrelation->Fill(fill);             
+           }
+         }             
+       }
+      }
+    }
+    if(unsetvtx) D0particle->UnsetOwnPrimaryVtx();
+  } // end loop on D*->Kpipi
+  
+  fHistEventsProcessed->Fill(0);
+  
+  PostData(1,fHistEventsProcessed) ;
+  PostData(2,fCFManager->GetParticleContainer()) ;
+  PostData(3,fCorrelation) ;
+}
+
+
+//___________________________________________________________________________
+void AliCFTaskForDStarAnalysis::Terminate(Option_t*)
+{
+  // The Terminate()   
+  AliAnalysisTaskSE::Terminate();
+  
+  // draw correlation matrix
+  AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
+  if(!cont) {
+    printf("CONTAINER NOT FOUND\n");
+    return;
+  } 
+}
+
+//___________________________________________________________________________
+void AliCFTaskForDStarAnalysis::UserCreateOutputObjects() {
+
+  //
+  // useroutput
+  //
+  
+  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
+  
+  //slot #1
+  OpenFile(1);
+  fHistEventsProcessed = new TH1I("CFDSchist0","",1,0,1) ;
+}
+
+//________________________________________________________________________________
+Bool_t AliCFTaskForDStarAnalysis::GetDStarMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC)const {
+  
+  // 
+  // fill the D* and D0 MC container
+  //
+  
+  Bool_t isDStar = kFALSE;
+  
+  if(TMath::Abs(mcPart->GetPdgCode())!=413) return isDStar;
+  
+  // 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 D* MC doesn't have correct daughters!!");
+    return isDStar;  
+  }
+  
+  if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
+    AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
+    return isDStar;  
+  }
+  
+  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
+  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
+  if (!mcPartDaughter0 || !mcPartDaughter1) {
+    AliWarning("D*: At least one Daughter Particle not found in tree, skipping"); 
+    return isDStar;  
+  }
+  
+  if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
+       TMath::Abs(mcPartDaughter1->GetPdgCode())==211) && 
+      !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
+       TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
+    AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
+    return isDStar;  
+  }
+  
+  Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
+  Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
+
+  // getting vertex from daughters
+  mcPartDaughter0->XvYvZv(vtx2daughter0);  
+  mcPartDaughter1->XvYvZv(vtx2daughter1);  
+  
+  // check if the secondary vertex is the same for both
+  if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
+    AliError("The D* daughters have different secondary vertex, skipping the track");
+    return isDStar;
+  }
+  
+  AliAODMCParticle* neutralDaugh = mcPartDaughter0;
+  
+  Double_t VectorD0[2] ={0.,0.};
+  
+  if (!EvaluateIfD0toKpi(neutralDaugh,mcArray,VectorD0)) {
+    AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
+    return isDStar;  
+  }
+  // get the pT of the daughters
+  
+  Double_t pTpi = 0.;
+  Double_t pTD0 = 0.;
+  
+  if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
+    pTpi = mcPartDaughter0->Pt();
+    pTD0 = mcPartDaughter1->Pt();
+  }
+  else {
+    pTpi = mcPartDaughter1->Pt();
+    pTD0 = mcPartDaughter0->Pt();
+  }
+  
+  vectorMC[0] = mcPart->Pt();
+  vectorMC[1] = mcPart->Y() ;
+  vectorMC[2] = 0;
+  vectorMC[3] = pTpi ;
+  vectorMC[4] = pTD0 ;
+  vectorMC[5] = 0;
+  vectorMC[6] = mcPart->Phi() ;
+  vectorMC[7] = VectorD0[0] ;
+  vectorMC[8] = VectorD0[1] ;
+  
+  isDStar = kTRUE;
+  
+  return isDStar;
+}
+//________________________________________________________________________________________________
+
+Bool_t AliCFTaskForDStarAnalysis::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray, Double_t* VectorD0)const{
+  
+  //
+  // chack wether D0 is decaing into kpi
+  //
+  
+  Bool_t isHadronic = kFALSE;
+  
+  Int_t daughterD00 = neutralDaugh->GetDaughter(0);
+  Int_t daughterD01 = neutralDaugh->GetDaughter(1);
+  
+  AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
+  if (daughterD00 == 0 || daughterD01 == 0) {
+    AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
+    return isHadronic;  
+  }
+  
+  if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
+    AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
+    return isHadronic;  
+  }
+  
+  AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
+  AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
+  if (!mcPartDaughterD00 || !mcPartDaughterD01) {
+    AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping"); 
+    return isHadronic;  
+  }
+  
+  if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
+       TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) && 
+      !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
+       TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
+    AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
+    return isHadronic;  
+  }
+
+  Double_t pTD0pi = 0;
+  Double_t pTD0K = 0;
+
+  if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
+    pTD0pi = mcPartDaughterD00->Pt();
+    pTD0K = mcPartDaughterD01->Pt();
+  }
+  else {
+    pTD0pi = mcPartDaughterD01->Pt();
+    pTD0K  = mcPartDaughterD00->Pt();
+  }
+  isHadronic = kTRUE;
+
+  VectorD0[0] = pTD0pi;
+  VectorD0[1] = pTD0K;
+
+  return isHadronic;
+
+}
diff --git a/PWG3/vertexingHF/AliCFTaskForDStarAnalysis.h b/PWG3/vertexingHF/AliCFTaskForDStarAnalysis.h
new file mode 100644 (file)
index 0000000..ef7f057
--- /dev/null
@@ -0,0 +1,90 @@
+#ifndef ALICFTASKFORDSTARANALYSIS_H
+#define ALICFTASKFORDSTARANALYSIS_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 D* corrections --  
+
+#include "AliAnalysisTaskSE.h"
+
+class TH1I;
+class TParticle ;
+class TFile ;
+class TClonesArray ;
+class AliCFManager;
+class AliAODRecoDecay;
+class AliAODRecoDecayHF2Prong;
+class AliAODMCParticle;
+class THnSparse;
+
+class AliCFTaskForDStarAnalysis : public AliAnalysisTaskSE {
+  public:
+
+  enum {
+    kStepGenerated       = 0,
+    kStepAcceptance      = 1,
+    kStepVertex          = 2,
+    kStepRefit           = 3,
+    kStepReconstructed   = 4,
+    kStepRecoAcceptance  = 5,
+    kStepRecoITSClusters = 6,
+    kStepRecoCuts        = 7
+  };
+
+  AliCFTaskForDStarAnalysis();
+  AliCFTaskForDStarAnalysis(const Char_t* name);
+  AliCFTaskForDStarAnalysis& operator= (const AliCFTaskForDStarAnalysis& c);
+  AliCFTaskForDStarAnalysis(const AliCFTaskForDStarAnalysis& c);
+  virtual ~AliCFTaskForDStarAnalysis();
+
+  // ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects
+  void     UserCreateOutputObjects();
+  void     UserExec(Option_t *option);
+  void     Terminate(Option_t *);
+
+ // UNFOLDING
+  void     SetCorrelationMatrix(THnSparse* h) {fCorrelation=h;}
+  void     SetAcceptanceUnf(Bool_t AcceptanceUnf) {fAcceptanceUnf = AcceptanceUnf;}
+  Bool_t   GetAcceptanceUnf() const {return fAcceptanceUnf;}
+
+  // CORRECTION FRAMEWORK
+  void           SetCFManager(AliCFManager* io) {fCFManager = io;}   // global correction manager
+  AliCFManager * GetCFManager()                 {return fCFManager;} // get corr manager
+
+  Bool_t   GetDStarMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC)const;
+  Bool_t   EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray, Double_t* VectorD0)const;
+  // for the D0 
+  void     SetMinITSClusters(Int_t minITSClusters) {fMinITSClusters = minITSClusters;}
+  Int_t    GetMinITSClusters() const {return fMinITSClusters;}
+  // for the soft pion
+  void     SetMinITSClustersSoft(Int_t minITSClustersSoft) {fMinITSClustersSoft = minITSClustersSoft;}
+  Int_t    GetMinITSClustersSoft() const {return fMinITSClustersSoft;}
+
+ protected:
+  AliCFManager   *fCFManager;   //  pointer to the CF manager
+  TH1I *fHistEventsProcessed;   //! simple histo for monitoring the number of events processed
+  THnSparse* fCorrelation;      //  response matrix for unfolding
+  Int_t fCountRecoDStarSel;     //  Reco particle found that satisfy cuts in D* selection
+  Int_t fEvents;                //  n. of events
+  Int_t fMinITSClusters;        //  min n. of ITS clusters for RecoDecay
+  Int_t fMinITSClustersSoft;    //  min n. of ITS clusters for RecoDecay soft pion
+  Bool_t fAcceptanceUnf;        //  flag for unfolding before or after cuts.
+  
+  ClassDef(AliCFTaskForDStarAnalysis,3); // class for D* corrections
+};
+
+#endif