New common interface to the CORRFW for all charm analyses (Davide, ChiaraZ)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 11 Sep 2010 01:16:38 +0000 (01:16 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 11 Sep 2010 01:16:38 +0000 (01:16 +0000)
PWG3/vertexingHF/AliCFTaskVertexingHF.cxx [new file with mode: 0644]
PWG3/vertexingHF/AliCFTaskVertexingHF.h [new file with mode: 0644]
PWG3/vertexingHF/AliCFVertexingHF.cxx [new file with mode: 0644]
PWG3/vertexingHF/AliCFVertexingHF.h [new file with mode: 0644]
PWG3/vertexingHF/AliCFVertexingHF2Prong.cxx [new file with mode: 0644]
PWG3/vertexingHF/AliCFVertexingHF2Prong.h [new file with mode: 0644]
PWG3/vertexingHF/macros/AddTaskCFVertexingHF.C [new file with mode: 0644]

diff --git a/PWG3/vertexingHF/AliCFTaskVertexingHF.cxx b/PWG3/vertexingHF/AliCFTaskVertexingHF.cxx
new file mode 100644 (file)
index 0000000..1eca2d4
--- /dev/null
@@ -0,0 +1,973 @@
+/**************************************************************************
+ * 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
+// 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl, 
+// Reco Acc + ITS Cl + PPR cuts
+// 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
+// dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
+//
+//-----------------------------------------------------------------------
+// Author : C. Zampolli, CERN
+//          D. Caffarri, Univ & INFN Padova  caffarri@pd.infn.it
+//-----------------------------------------------------------------------
+//-----------------------------------------------------------------------
+// Base class for HF Unfolding (pt and eta)
+// correlation matrix filled at Acceptance and PPR level
+// Author: A.Grelli ,  Utrecht - agrelli@uu.nl
+//----------------------------------------------------------------------- 
+#include <TCanvas.h>
+#include <TParticle.h>
+#include <TDatabasePDG.h>
+#include <TH1I.h>
+#include <TStyle.h>
+#include <TFile.h>
+
+#include "AliCFTaskVertexingHF.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 "AliAODRecoDecayHF2Prong.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "AliESDtrack.h"
+#include "TChain.h"
+#include "THnSparse.h"
+#include "TH2D.h"
+#include "AliESDtrackCuts.h"
+#include "AliRDHFCutsD0toKpi.h"
+#include "AliCFVertexingHF2Prong.h"
+#include "AliCFVertexingHF.h"
+
+//__________________________________________________________________________
+AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
+AliAnalysisTaskSE(),
+fPDG(0),
+fCFManager(0x0),
+fHistEventsProcessed(0x0),
+fCorrelation(0x0),
+fCountMC(0),
+fCountAcc(0),
+fCountVertex(0),
+fCountRefit(0),
+fCountReco(0),
+fCountRecoAcc(0),
+fCountRecoITSClusters(0),
+fCountRecoPPR(0),
+fCountRecoPID(0),
+fEvents(0),
+fFillFromGenerated(kFALSE),
+fOriginDselection(0),
+fAcceptanceUnf(kTRUE),
+fCuts(0)
+{
+       //
+       //Default ctor
+       //
+}
+//___________________________________________________________________________
+AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCutsD0toKpi* cuts) :
+AliAnalysisTaskSE(name),
+fPDG(0),
+fCFManager(0x0),
+fHistEventsProcessed(0x0),
+fCorrelation(0x0),
+fCountMC(0),
+fCountAcc(0),
+fCountVertex(0),
+fCountRefit(0),
+fCountReco(0),
+fCountRecoAcc(0),
+fCountRecoITSClusters(0),
+fCountRecoPPR(0),
+fCountRecoPID(0),
+fEvents(0),
+fFillFromGenerated(kFALSE),
+fOriginDselection(0),
+fAcceptanceUnf(kTRUE),
+fCuts(cuts)
+{
+       //
+       // Constructor. Initialization of Inputs and Outputs
+       //
+       Info("AliCFTaskVertexingHF","Calling Constructor");
+       /*
+        DefineInput(0) and DefineOutput(0)
+        are taken care of by AliAnalysisTaskSE constructor
+        */
+       DefineOutput(1,TH1I::Class());
+       DefineOutput(2,AliCFContainer::Class());
+       DefineOutput(3,THnSparseD::Class());
+       
+       fCuts->PrintAll();
+}
+
+//___________________________________________________________________________
+AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c) 
+{
+       //
+       // Assignment operator
+       //
+       if (this!=&c) {
+               AliAnalysisTaskSE::operator=(c) ;
+               fPDG      = c.fPDG;
+               fCFManager  = c.fCFManager;
+               fHistEventsProcessed = c.fHistEventsProcessed;
+               fCuts = c.fCuts;
+       }
+       return *this;
+}
+
+//___________________________________________________________________________
+AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
+AliAnalysisTaskSE(c),
+fPDG(c.fPDG),
+fCFManager(c.fCFManager),
+fHistEventsProcessed(c.fHistEventsProcessed),
+fCorrelation(c.fCorrelation),
+fCountMC(c.fCountMC),
+fCountAcc(c.fCountAcc),
+fCountVertex(c.fCountVertex),
+fCountRefit(c.fCountRefit),
+fCountReco(c.fCountReco),
+fCountRecoAcc(c.fCountRecoAcc),
+fCountRecoITSClusters(c.fCountRecoITSClusters),
+fCountRecoPPR(c.fCountRecoPPR),
+fCountRecoPID(c.fCountRecoPID),
+fEvents(c.fEvents),
+fFillFromGenerated(c.fFillFromGenerated),
+fOriginDselection(c.fOriginDselection),
+fAcceptanceUnf(c.fAcceptanceUnf),
+fCuts(c.fCuts)
+{
+       //
+       // Copy Constructor
+       //
+}
+
+//___________________________________________________________________________
+AliCFTaskVertexingHF::~AliCFTaskVertexingHF() {
+       //
+       //destructor
+       //
+       if (fCFManager)           delete fCFManager ;
+       if (fHistEventsProcessed) delete fHistEventsProcessed ;
+       if (fCorrelation)                 delete fCorrelation ;
+       if (fCuts)                delete fCuts;
+}
+
+//_________________________________________________
+void AliCFTaskVertexingHF::UserExec(Option_t *)
+{
+       
+       //
+       // Main loop function
+       //
+       
+       PostData(1,fHistEventsProcessed) ;
+       PostData(2,fCFManager->GetParticleContainer()) ;
+       PostData(3,fCorrelation) ;
+       
+       AliESDtrackCuts* trackCuts = fCuts->GetTrackCuts();
+       
+       if (fFillFromGenerated){
+               AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
+       }
+       
+       if (!fInputEvent) {
+               Error("UserExec","NO EVENT FOUND!");
+               return;
+       }
+       
+       AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
+       
+       TClonesArray *arrayD0toKpi=0;
+       //TClonesArray *arrayBranch=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();
+           arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+         }
+       } else {
+         arrayD0toKpi=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
+       }
+       
+       AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
+       if (!aodVtx) return;
+
+       if (!arrayD0toKpi) {
+         AliError("Could not find array of HF vertices");
+         return;
+       }
+       
+       fEvents++;
+       if (fEvents%10000 == 0) AliDebug(2,Form("Event %d",fEvents));
+       
+       fCFManager->SetRecEventInfo(aodEvent);
+       fCFManager->SetMCEventInfo(aodEvent);
+       
+       
+       //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
+       Int_t nVar = 13;
+       
+       Double_t* containerInput = new Double_t[nVar];
+       Double_t* containerInputMC = new Double_t[nVar]; 
+       
+       //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;
+       }
+       Int_t icountMC = 0;
+       Int_t icountAcc = 0;
+       Int_t icountReco = 0;
+       Int_t icountVertex = 0;
+       Int_t icountRefit = 0;
+       Int_t icountRecoAcc = 0;
+       Int_t icountRecoITSClusters = 0;
+       Int_t icountRecoPPR = 0;
+       Int_t icountRecoPID = 0;
+       
+       Int_t cquarks = 0;
+
+       // These flgs will be setted i the config macro
+       UShort_t originDselection = 0;
+       
+       
+       AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+       if (!mcHeader) {
+               AliError("Could not find MC Header in AOD");
+               return;
+       }
+
+
+       AliCFVertexingHF2Prong *cfVtxHF = new AliCFVertexingHF2Prong(mcArray, originDselection);
+
+       Double_t zPrimVertex = aodVtx ->GetZ();
+       Double_t zMCVertex = mcHeader->GetVtxZ();
+
+       //General settings: vertex, feed down and fill reco container with generated values.                    
+       cfVtxHF->SetRecoPrimVertex(zPrimVertex);
+       cfVtxHF->SetMCPrimaryVertex(zMCVertex);
+       cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
+
+       //      cfVtxHF->SetKeepD0fromBOnly(fKeepD0fromBOnly);
+               //      cfVtxHF->SetKeepD0fromB(fKeepD0fromB);
+
+       for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { 
+
+         AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
+
+         // check the MC-level cuts, must be a D0
+         if (!fCFManager->CheckParticleCuts(0, mcPart)) continue;  // 0 stands for MC level
+         
+         cfVtxHF->SetMCCandidateParam(iPart);
+         cfVtxHF->SetNVar(nVar);
+         //counting c quarks
+         cquarks += cfVtxHF->MCcquarkCounting(mcPart);
+        
+         //check the candiate family at MV level
+         if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) continue;
+
+         //Fill the MC container
+         Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
+         
+         if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
+         if (TMath::Abs(containerInputMC[1]) < 0.5) {
+           fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc);
+         }
+           
+         if (mcContainerFilled) {
+           
+           //MC 
+           fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated);
+           icountMC++;
+           printf("MC cointainer filled \n");
+           
+             // MC in acceptance
+             // check the MC-Acceptance level cuts
+             // since standard CF functions are not applicable, using Kine Cuts on daughters
+             Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
+             if (mcAccepStep){ 
+                         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
+                         printf("MC acceptance cut passed\n");
+                         icountAcc++;
+               
+                         //MC Vertex step
+                         if (fCuts->IsEventSelected(aodEvent)){
+                                 // filling the container if the vertex is ok
+                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
+                                 printf("Vertex cut passed and container filled\n");
+                                 icountVertex++;
+                 
+                                 //mc Refit requirement        
+                                 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, trackCuts);
+                                 if (mcRefitStep){
+                                         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
+                                         printf("MC Refit cut passed and container filled\n");
+                                         icountRefit++;
+                                 }
+                                 else{
+                                         AliDebug(3,"MC Refit cut not passed\n");
+                                         continue;
+                                 }
+                 
+                         }
+                         else{
+                                 AliDebug (3, "MC vertex step not passed\n");
+                                 continue;
+                         }
+                 }
+                 else{
+                         AliDebug (3, "MC in acceptance step not passed\n");
+                         continue;
+                 }
+               
+         }
+         else {
+                 AliDebug (3, "MC container not filled\n");
+               }
+       }
+               
+       if (cquarks<2) AliDebug(2,Form("Eveith %d c-quarks", cquarks));
+       AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
+       AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
+       AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
+       AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
+       
+       // Now go to rec level
+       fCountMC += icountMC;
+       fCountAcc += icountAcc;
+       fCountVertex+= icountVertex;
+       fCountRefit+= icountRefit;
+       
+       AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
+       
+
+       for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
+         
+         AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
+               Bool_t unsetvtx=kFALSE;
+               if(!d0tokpi->GetOwnPrimaryVtx()) {
+                       d0tokpi->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
+                       unsetvtx=kTRUE;
+               }
+         
+               if (d0tokpi->GetPrimaryVtx()) printf("candidate has prim vtx \n");       
+               if (aodEvent->GetPrimaryVertex()) printf ("aod event vertex\n");
+         
+               Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)d0tokpi);
+               if (!signAssociation){
+                       d0tokpi = 0x0;
+                       continue;
+               }
+       
+               Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
+               if (recoContFilled){
+
+                       if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;         
+
+                       //Reco Step
+                       Bool_t recoStep = cfVtxHF->RecoStep();
+                       Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
+           
+                       if (recoStep && recoContFilled && vtxCheck){
+                               fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
+                               icountReco++;
+                               printf("Reco step  passed and container filled\n");
+           
+           
+                               //Reco in the acceptance -- take care of UNFOLDING!!!!
+                               Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(trackCuts);
+                               if (recoAcceptanceStep) {
+                                       fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
+                                       icountRecoAcc++; 
+                                       printf("Reco acceptance cut passed and container filled\n");
+               
+                                       if(fAcceptanceUnf){
+                                               Double_t fill[4]; //fill response matrix
+                                               Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
+                                               if (bUnfolding) fCorrelation->Fill(fill);
+                                               }
+               
+                                       //Number of ITS cluster requirements    
+                                       Int_t recoITSnCluster = fCuts->IsSelected(d0tokpi, AliRDHFCuts::kTracks);
+                                       if (recoITSnCluster){
+                                               fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
+                                               icountRecoITSClusters++;   
+                                               printf("Reco n ITS cluster cut passed and container filled\n");
+           
+
+                                               Int_t recoAnalysisCuts = fCuts->IsSelected(d0tokpi, AliRDHFCuts::kCandidate);
+                                               if(recoAnalysisCuts){
+                                                 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR);
+                                                 icountRecoPPR++;
+                                                 printf ("Reco Analysis cuts passed and container filled \n");
+                                                 
+                                                 
+                                                 //pid selection       
+                                                 Int_t recoPidSelection = fCuts->IsSelected(d0tokpi, AliRDHFCuts::kPID);
+                                                 if (recoPidSelection > 0){
+                                                   fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID);
+                                                   icountRecoPID++;
+                                                   printf ("Reco PID cuts passed and container filled \n");
+                                                   
+                                                   
+                                                   if(!fAcceptanceUnf){
+                                                     Double_t fill[4]; //fill response matrix
+                                                     Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
+                                                     if (bUnfolding) fCorrelation->Fill(fill);
+                                                   }
+                                                 }
+                                                 else {
+                                                   AliDebug(3, "Analysis Cuts step not passed \n");
+                                                   continue;
+                                                 }
+                                               }
+                                               else {
+                                                 AliDebug(3, "PID selection not passed \n");
+                                                 continue;
+                                               }
+                                       }
+                                       else{
+                                         AliDebug(3, "Number of ITS cluster step not passed\n");
+                                         continue;
+                                       }
+                               }
+                               else{
+                                 AliDebug(3, "Reco acceptance step not passed\n");
+                                 continue;
+                               }
+                       }
+                       else {
+                         AliDebug(3, "Reco step not passed\n");
+                         continue;
+                       }
+               }
+               
+               if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
+       } // end loop on D0->Kpi
+       
+       AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
+       
+       fCountReco+= icountReco;
+               fCountRecoAcc+= icountRecoAcc;
+       fCountRecoITSClusters+= icountRecoITSClusters;
+       fCountRecoPPR+= icountRecoPPR;
+       fCountRecoPID+= icountRecoPID;
+       
+       fHistEventsProcessed->Fill(0);
+}
+
+//___________________________________________________________________________
+void AliCFTaskVertexingHF::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.
+       
+       AliAnalysisTaskSE::Terminate();
+       
+       AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
+       AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
+       AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
+       AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fEvents));
+       AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
+       AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
+       AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fEvents));
+       AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
+       
+       // draw some example plots....
+       
+       AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
+       if(!cont) {
+               printf("CONTAINER NOT FOUND\n");
+               return;
+       }
+       // projecting the containers to obtain histograms
+       // first argument = variable, second argument = step
+       
+       // MC-level
+       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* h60 =   cont->ShowProjection(6,0) ;   // dca
+       TH1D* h70 =   cont->ShowProjection(7,0) ;   // d0pi
+       TH1D* h80 =   cont->ShowProjection(8,0) ;   // d0K
+       TH1D* h90 =   cont->ShowProjection(9,0) ;   // d0xd0
+       TH1D* h100 =   cont->ShowProjection(10,0) ;   // cosPointingAngle
+       TH1D* h110 =   cont->ShowProjection(11,0) ;   // phi
+       
+       // MC-Acceptance level
+       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
+       TH1D* h61 =   cont->ShowProjection(6,1) ;   // dca
+       TH1D* h71 =   cont->ShowProjection(7,1) ;   // d0pi
+       TH1D* h81 =   cont->ShowProjection(8,1) ;   // d0K
+       TH1D* h91 =   cont->ShowProjection(9,1) ;   // d0xd0
+       TH1D* h101 =   cont->ShowProjection(10,1) ;   // cosPointingAngle
+       TH1D* h111 =   cont->ShowProjection(11,1) ;   // phi
+       
+       // Reco-level
+       TH1D* h04 =   cont->ShowProjection(0,4) ;   // pt
+       TH1D* h14 =   cont->ShowProjection(1,4) ;   // rapidity
+       TH1D* h24 =   cont->ShowProjection(2,4) ;   // cosThetaStar
+       TH1D* h34 =   cont->ShowProjection(3,4) ;   // pTpi
+       TH1D* h44 =   cont->ShowProjection(4,4) ;   // pTK
+       TH1D* h54 =   cont->ShowProjection(5,4) ;   // cT
+       TH1D* h64 =   cont->ShowProjection(6,4) ;   // dca
+       TH1D* h74 =   cont->ShowProjection(7,4) ;   // d0pi
+       TH1D* h84 =   cont->ShowProjection(8,4) ;   // d0K
+       TH1D* h94 =   cont->ShowProjection(9,4) ;   // d0xd0
+       TH1D* h104 =   cont->ShowProjection(10,4) ;   // cosPointingAngle
+       TH1D* h114 =   cont->ShowProjection(11,4) ;   // phi
+       
+       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)");
+       h60->SetTitle("dca (#mum)");
+       h70->SetTitle("d0_pi (#mum)");
+       h80->SetTitle("d0_K (#mum)");
+       h90->SetTitle("d0xd0 (#mum^2)");
+       h100->SetTitle("cosPointingAngle");
+       h100->SetTitle("phi (rad)");
+       
+       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)");
+       h60->GetXaxis()->SetTitle("dca (#mum)");
+       h70->GetXaxis()->SetTitle("d0_pi (#mum)");
+       h80->GetXaxis()->SetTitle("d0_K (#mum)");
+       h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
+       h100->GetXaxis()->SetTitle("cosPointingAngle");
+       h110->GetXaxis()->SetTitle("phi (rad)");
+       
+       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)");
+       h61->SetTitle("dca (#mum)");
+       h71->SetTitle("d0_pi (#mum)");
+       h81->SetTitle("d0_K (#mum)");
+       h91->SetTitle("d0xd0 (#mum^2)");
+       h101->SetTitle("cosPointingAngle");
+       h111->GetXaxis()->SetTitle("phi (rad)");
+       
+       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)");
+       h61->GetXaxis()->SetTitle("dca (#mum)");
+       h71->GetXaxis()->SetTitle("d0_pi (#mum)");
+       h81->GetXaxis()->SetTitle("d0_K (#mum)");
+       h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
+       h101->GetXaxis()->SetTitle("cosPointingAngle");
+       h111->GetXaxis()->SetTitle("phi (rad)");
+       
+       h04->SetTitle("pT_D0 (GeV/c)");
+       h14->SetTitle("rapidity");
+       h24->SetTitle("cosThetaStar");
+       h34->SetTitle("pT_pi (GeV/c)");
+       h44->SetTitle("pT_K (Gev/c)");
+       h54->SetTitle("cT (#mum)");
+       h64->SetTitle("dca (#mum)");
+       h74->SetTitle("d0_pi (#mum)");
+       h84->SetTitle("d0_K (#mum)");
+       h94->SetTitle("d0xd0 (#mum^2)");
+       h104->SetTitle("cosPointingAngle");
+       h114->GetXaxis()->SetTitle("phi (rad)");
+       
+       h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
+       h14->GetXaxis()->SetTitle("rapidity");
+       h24->GetXaxis()->SetTitle("cosThetaStar");
+       h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
+       h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
+       h54->GetXaxis()->SetTitle("cT (#mum)");
+       h64->GetXaxis()->SetTitle("dca (#mum)");
+       h74->GetXaxis()->SetTitle("d0_pi (#mum)");
+       h84->GetXaxis()->SetTitle("d0_K (#mum)");
+       h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
+       h104->GetXaxis()->SetTitle("cosPointingAngle");
+       h114->GetXaxis()->SetTitle("phi (rad)");
+       
+       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();
+       Double_t max6 = h60->GetMaximum();
+       Double_t max7 = h70->GetMaximum();
+       Double_t max8 = h80->GetMaximum();
+       Double_t max9 = h90->GetMaximum();
+       Double_t max10 = h100->GetMaximum();
+       Double_t max11 = h110->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);
+       h60->GetYaxis()->SetRangeUser(0,max6*1.2);
+       h70->GetYaxis()->SetRangeUser(0,max7*1.2);
+       h80->GetYaxis()->SetRangeUser(0,max8*1.2);
+       h90->GetYaxis()->SetRangeUser(0,max9*1.2);
+       h100->GetYaxis()->SetRangeUser(0,max10*1.2);
+       h110->GetYaxis()->SetRangeUser(0,max11*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);
+       h61->GetYaxis()->SetRangeUser(0,max6*1.2);
+       h71->GetYaxis()->SetRangeUser(0,max7*1.2);
+       h81->GetYaxis()->SetRangeUser(0,max8*1.2);
+       h91->GetYaxis()->SetRangeUser(0,max9*1.2);
+       h101->GetYaxis()->SetRangeUser(0,max10*1.2);
+       h111->GetYaxis()->SetRangeUser(0,max11*1.2);
+       
+       h00->SetMarkerStyle(20);
+       h10->SetMarkerStyle(24);
+       h20->SetMarkerStyle(21);
+       h30->SetMarkerStyle(25);
+       h40->SetMarkerStyle(27);
+       h50->SetMarkerStyle(28);
+       h60->SetMarkerStyle(20);
+       h70->SetMarkerStyle(24);
+       h80->SetMarkerStyle(21);
+       h90->SetMarkerStyle(25);
+       h100->SetMarkerStyle(27);
+       h110->SetMarkerStyle(28);
+       
+       h00->SetMarkerColor(2);
+       h10->SetMarkerColor(2);
+       h20->SetMarkerColor(2);
+       h30->SetMarkerColor(2);
+       h40->SetMarkerColor(2);
+       h50->SetMarkerColor(2);
+       h60->SetMarkerColor(2);
+       h70->SetMarkerColor(2);
+       h80->SetMarkerColor(2);
+       h90->SetMarkerColor(2);
+       h100->SetMarkerColor(2);
+       h110->SetMarkerColor(2);
+       
+       h01->SetMarkerStyle(20) ;
+       h11->SetMarkerStyle(24) ;
+       h21->SetMarkerStyle(21) ;
+       h31->SetMarkerStyle(25) ;
+       h41->SetMarkerStyle(27) ;
+       h51->SetMarkerStyle(28) ;
+       h61->SetMarkerStyle(20);
+       h71->SetMarkerStyle(24);
+       h81->SetMarkerStyle(21);
+       h91->SetMarkerStyle(25);
+       h101->SetMarkerStyle(27);
+       h111->SetMarkerStyle(28);
+       
+       h01->SetMarkerColor(8);
+       h11->SetMarkerColor(8);
+       h21->SetMarkerColor(8);
+       h31->SetMarkerColor(8);
+       h41->SetMarkerColor(8);
+       h51->SetMarkerColor(8);
+       h61->SetMarkerColor(8);
+       h71->SetMarkerColor(8);
+       h81->SetMarkerColor(8);
+       h91->SetMarkerColor(8);
+       h101->SetMarkerColor(8);
+       h111->SetMarkerColor(8);
+       
+       h04->SetMarkerStyle(20) ;
+       h14->SetMarkerStyle(24) ;
+       h24->SetMarkerStyle(21) ;
+       h34->SetMarkerStyle(25) ;
+       h44->SetMarkerStyle(27) ;
+       h54->SetMarkerStyle(28) ;
+       h64->SetMarkerStyle(20);
+       h74->SetMarkerStyle(24);
+       h84->SetMarkerStyle(21);
+       h94->SetMarkerStyle(25);
+       h104->SetMarkerStyle(27);
+       h114->SetMarkerStyle(28);
+       
+       h04->SetMarkerColor(4);
+       h14->SetMarkerColor(4);
+       h24->SetMarkerColor(4);
+       h34->SetMarkerColor(4);
+       h44->SetMarkerColor(4);
+       h54->SetMarkerColor(4);
+       h64->SetMarkerColor(4);
+       h74->SetMarkerColor(4);
+       h84->SetMarkerColor(4);
+       h94->SetMarkerColor(4);
+       h104->SetMarkerColor(4);
+       h114->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("c1New","pT, rapidiy, cosThetaStar",1100,1600);
+       c1->Divide(3,3);
+       
+       c1->cd(1);
+       h00->Draw("p");
+       c1->cd(1);
+       c1->cd(2);
+       h01->Draw("p");
+       c1->cd(2);
+       c1->cd(3);
+       h04->Draw("p");
+       c1->cd(3);
+       c1->cd(4);
+       h10->Draw("p");
+       c1->cd(4);
+       c1->cd(5);
+       h11->Draw("p");
+       c1->cd(5);
+       c1->cd(6);
+       h14->Draw("p");
+       c1->cd(6);
+       c1->cd(7);
+       h20->Draw("p");
+       c1->cd(7);
+       c1->cd(8);
+       h21->Draw("p");
+       c1->cd(8);
+       c1->cd(9);
+       h24->Draw("p");
+       c1->cd(9);
+       c1->cd();
+       
+       TCanvas * c2 =new TCanvas("c2New","pTpi, pTK, cT",1100,1600);
+       c2->Divide(3,3);
+       c2->cd(1);
+       h30->Draw("p");
+       c2->cd(1);
+       c2->cd(2);
+       h31->Draw("p");
+       c2->cd(2);
+       c2->cd(3);
+       h34->Draw("p");
+       c2->cd(3);
+       c2->cd(4);
+       h40->Draw("p");
+       c2->cd(4);
+       c2->cd(5);
+       h41->Draw("p");
+       c2->cd(5);
+       c2->cd(6);
+       h44->Draw("p");
+       c2->cd(6);
+       c2->cd(7);
+       h50->Draw("p");
+       c2->cd(7);
+       c2->cd(8);
+       h51->Draw("p");
+       c2->cd(8);
+       c2->cd(9);
+       h54->Draw("p");
+       c2->cd(9);
+       c2->cd();
+       
+       TCanvas * c3 =new TCanvas("c3New","dca, d0pi, d0K",1100,1600);
+       c3->Divide(3,3);
+       c3->cd(1);
+       h60->Draw("p");
+       c3->cd(1);
+       c3->cd(2);
+       h61->Draw("p");
+       c3->cd(2);
+       c3->cd(3);
+       h64->Draw("p");
+       c3->cd(3);
+       c3->cd(4);
+       h70->Draw("p");
+       c3->cd(4);
+       c3->cd(5);
+       h71->Draw("p");
+       c3->cd(5);
+       c3->cd(6);
+       h74->Draw("p");
+       c3->cd(6);
+       c3->cd(7);
+       h80->Draw("p");
+       c3->cd(7);
+       c3->cd(8);
+       h81->Draw("p");
+       c3->cd(8);
+       c3->cd(9);
+       h84->Draw("p");
+       c3->cd(9);
+       c3->cd();
+       
+       TCanvas * c4 =new TCanvas("c4New","d0xd0, cosPointingAngle, phi",1100,1600);
+       c4->Divide(3,3);
+       c4->cd(1);
+       h90->Draw("p");
+       c4->cd(1);
+       c4->cd(2);
+       h91->Draw("p");
+       c4->cd(2);
+       c4->cd(3);
+       h94->Draw("p");
+       c4->cd(3);
+       c4->cd(4);
+       h100->Draw("p");
+       c4->cd(4);
+       c4->cd(5);
+       h101->Draw("p");
+       c4->cd(5);
+       c4->cd(6);
+       h104->Draw("p");
+       c4->cd(6);
+       c4->cd(7);
+       h110->Draw("p");
+       c4->cd(7);
+       c4->cd(8);
+       h111->Draw("p");
+       c4->cd(8);
+       c4->cd(9);
+       h114->Draw("p");
+       c4->cd(9);
+       c4->cd();
+       
+       THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
+       
+       TH2D* corr1 =hcorr->Projection(0,2);
+       TH2D* corr2 = hcorr->Projection(1,3);
+       
+       TCanvas * c7 =new TCanvas("c7New","",800,400);
+       c7->Divide(2,1);
+       c7->cd(1);
+       corr1->Draw("text");
+       c7->cd(2);
+       corr2->Draw("text");
+       
+       
+       TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
+       
+       corr1->Write();
+       corr2->Write();
+       h00->Write("pT_D0_step0");
+       h10->Write("rapidity_step0");
+       h20->Write("cosThetaStar_step0");
+       h30->Write("pT_pi_step0");
+       h40->Write("pT_K_step0");
+       h50->Write("cT_step0");
+       h60->Write("dca_step0");
+       h70->Write("d0_pi_step0");
+       h80->Write("d0_K_step0");
+       h90->Write("d0xd0_step0");
+       h100->Write("cosPointingAngle_step0");
+       h110->Write("phi_step0");
+       
+       h01->Write("pT_D0_step1");
+       h11->Write("rapidity_step1");
+       h21->Write("cosThetaStar_step1");
+       h31->Write("pT_pi_step1");
+       h41->Write("pT_K_step1");
+       h51->Write("cT_step1");
+       h61->Write("dca_step1");
+       h71->Write("d0_pi_step1");
+       h81->Write("d0_K_step1");
+       h91->Write("d0xd0_step1");
+       h101->Write("cosPointingAngle_step1");
+       h111->Write("phi_step1");
+       
+       h04->Write("pT_D0_step2");
+       h14->Write("rapidity_step2");
+       h24->Write("cosThetaStar_step2");
+       h34->Write("pT_pi_step2");
+       h44->Write("pT_K_step2");
+       h54->Write("cT_step2");
+       h64->Write("dca_step2");
+       h74->Write("d0_pi_step2");
+       h80->Write("d0_K_step2");
+       h94->Write("d0xd0_step2");
+       h104->Write("cosPointingAngle_step2");
+       h114->Write("phi_step2");
+       
+       file_projection->Close();
+       
+    
+       
+       /*
+        c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
+        c2->SaveAs("Plots/pTpi_pTK_cT.eps");
+        c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
+        c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
+        
+        c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
+        c2->SaveAs("Plots/pTpi_pTK_cT.gif");
+        c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
+        c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
+        */     
+}
+
+//___________________________________________________________________________
+void AliCFTaskVertexingHF::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("CFHFchist0","",1,0,1) ;
+}
+
diff --git a/PWG3/vertexingHF/AliCFTaskVertexingHF.h b/PWG3/vertexingHF/AliCFTaskVertexingHF.h
new file mode 100644 (file)
index 0000000..d3d0f82
--- /dev/null
@@ -0,0 +1,114 @@
+#ifndef ALICFTASKVERTEXINGHF_H
+#define ALICFTASKVERTEXINGHF_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 and step 
+// Author : C. Zampolli, CERN
+//                     D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
+// Base class for HF Unfolding - agrelli@uu.nl
+//-----------------------------------------------------------------------
+
+
+#include "AliAnalysisTaskSE.h"
+#include "AliCFVertexingHF2Prong.h"
+#include "AliCFVertexingHF.h"
+
+class TH1I;
+class TParticle ;
+class TFile ;
+class TClonesArray ;
+class AliCFManager;
+class AliAODRecoDecay;
+class AliAODRecoDecayHF2Prong;
+class AliAODMCParticle;
+class THnSparse;
+class AliRDHFCutsD0toKpi;
+class AliCFVertexingHF2Prong;
+
+class AliCFTaskVertexingHF: public AliAnalysisTaskSE {
+public:
+       
+       enum {
+               kStepGeneratedLimAcc = 0,
+               kStepGenerated       = 1,
+               kStepAcceptance      = 2,
+               kStepVertex          = 3,
+               kStepRefit           = 4,
+               kStepReconstructed   = 5,
+               kStepRecoAcceptance  = 6,
+               kStepRecoITSClusters = 7,
+               kStepRecoPPR         = 8,
+               kStepRecoPID         = 9
+       };
+       
+       AliCFTaskVertexingHF();
+       AliCFTaskVertexingHF(const Char_t* name, AliRDHFCutsD0toKpi* cuts);
+       AliCFTaskVertexingHF& operator= (const AliCFTaskVertexingHF& c);
+       AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c);
+       virtual ~AliCFTaskVertexingHF();
+       
+       // 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 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
+
+       //      Setters for teh config macro
+       void   SetFillFromGenerated(Bool_t flag) {fFillFromGenerated = flag;}
+       Bool_t GetFillFromGenerated() const {return fFillFromGenerated;}
+
+
+       //void   SetDselection(UShort_t originDselection);
+        
+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
+       THnSparse* fCorrelation;      //  response matrix for unfolding
+       Int_t fCountMC;               //  MC particle found
+       Int_t fCountAcc;              //  MC particle found that satisfy acceptance cuts
+       Int_t fCountVertex;       //  Reco particle found that satisfy vertex constrained
+       Int_t fCountRefit;        //  Reco particle found that satisfy kTPCrefit and kITSrefit
+       Int_t fCountReco;             //  Reco particle found that satisfy cuts
+       Int_t fCountRecoAcc;          //  Reco particle found that satisfy cuts in requested acceptance
+       Int_t fCountRecoITSClusters;  //  Reco particle found that satisfy cuts in n. of ITS clusters
+       Int_t fCountRecoPPR;          //  Reco particle found that satisfy cuts in PPR
+       Int_t fCountRecoPID;          //Reco PID step 
+       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
+       UShort_t fOriginDselection;      // flag to select D0 origins. 0 Only from charm 1 only from beauty 2 both from charm and beauty
+       Bool_t fAcceptanceUnf;        //  flag for unfolding before or after cuts.
+       AliRDHFCutsD0toKpi* fCuts;    // cuts
+
+       
+       ClassDef(AliCFTaskVertexingHF,1); // class for HF corrections as a function of many variables
+};
+
+#endif
diff --git a/PWG3/vertexingHF/AliCFVertexingHF.cxx b/PWG3/vertexingHF/AliCFVertexingHF.cxx
new file mode 100644 (file)
index 0000000..f54ef64
--- /dev/null
@@ -0,0 +1,546 @@
+/**************************************************************************
+ * 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 and step 
+// Author : C. Zampolli, CERN
+// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
+// Base class for HF Unfolding - agrelli@uu.nl
+//-----------------------------------------------------------------------
+
+#include "TParticle.h"
+#include "TClonesArray.h"
+#include "AliAODMCParticle.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODRecoDecayHF3Prong.h"
+#include "AliAODRecoDecayHF4Prong.h"
+#include "AliAODMCHeader.h"
+#include "AliAODEvent.h"
+#include "AliLog.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDtrack.h"
+
+#include "AliCFVertexingHF.h"
+
+//___________________________________________________________
+AliCFVertexingHF::AliCFVertexingHF() :
+       fmcArray(0x0),
+       fRecoCandidate(0),
+       fmcPartCandidate(0x0),
+       fNDaughters(0),
+       fNVar(0),
+       fzPrimVertex(0),
+       fzMCVertex(0),
+       fFillFromGenerated(0),
+       fOriginDselection(0),
+       fKeepDfromB(kFALSE),
+       fKeepDfromBOnly(kFALSE),
+       fmcLabel(0),
+       fProngs(-1)
+{
+       return;
+}
+
+
+
+//_____________________________________________________
+AliCFVertexingHF::AliCFVertexingHF(TClonesArray *mcArray, UShort_t originDselection) :
+       fmcArray(mcArray),
+       fRecoCandidate(0),
+       fmcPartCandidate(0x0),
+       fNDaughters(0),
+       fNVar(0),
+       fzPrimVertex(0),
+       fzMCVertex(0),
+       fFillFromGenerated(0),
+       fOriginDselection(0),
+       fKeepDfromB(kFALSE),
+       fKeepDfromBOnly(kFALSE),
+       fmcLabel(0),
+       fProngs(-1)
+{
+
+  
+  SetDselection(originDselection);
+
+
+  return;
+}
+
+
+
+//_______________________________________________________
+AliCFVertexingHF::~AliCFVertexingHF()
+{
+       if (fmcArray) delete fmcArray;
+       if (fRecoCandidate) delete fRecoCandidate;
+       if (fmcPartCandidate) delete fmcPartCandidate;
+}
+
+
+//_____________________________________________________
+AliCFVertexingHF& AliCFVertexingHF::operator=(const AliCFVertexingHF& c){
+       
+       if (this!= &c){
+               TObject::operator=(c);
+               fmcArray = c.fmcArray;
+               fRecoCandidate = c.fRecoCandidate;
+               fmcPartCandidate = c.fmcPartCandidate;
+               fNDaughters = c.fNDaughters;
+               fNVar = c.fNVar;
+               fzPrimVertex = c.fzPrimVertex;
+               fzMCVertex = c.fzMCVertex;
+               fFillFromGenerated = c.fFillFromGenerated;
+               fOriginDselection = c.fOriginDselection;
+               fKeepDfromB = c.fKeepDfromB;
+               fKeepDfromBOnly = c.fKeepDfromBOnly;
+               fmcLabel = c.fmcLabel;
+               
+       }
+       
+       return *this;
+}
+
+//____________________________________________________
+AliCFVertexingHF::AliCFVertexingHF(const AliCFVertexingHF &c) :
+        TObject(c),
+       fmcArray(c.fmcArray),
+       fRecoCandidate(c.fRecoCandidate),
+       fmcPartCandidate(c.fmcPartCandidate),
+       fNDaughters(c.fNDaughters),
+       fNVar(c.fNVar),
+       fzPrimVertex(c.fzPrimVertex),
+       fzMCVertex(c.fzMCVertex),
+       fFillFromGenerated(c.fFillFromGenerated),
+       fOriginDselection (c.fOriginDselection),
+       fKeepDfromB (c.fKeepDfromB),
+       fKeepDfromBOnly (c.fKeepDfromBOnly),
+       fmcLabel(c.fmcLabel),
+       fProngs(c.fProngs)
+{
+   
+
+}
+
+//___________________________________________________________
+void AliCFVertexingHF::SetDselection(UShort_t originDselection){
+
+ fOriginDselection = originDselection;
+ if (fOriginDselection == 0) {
+   fKeepDfromB = kFALSE;
+   fKeepDfromBOnly = kFALSE;
+ }
+ if (fOriginDselection == 1) {
+   fKeepDfromB = kTRUE;
+   fKeepDfromBOnly = kTRUE;
+ }
+ if (fOriginDselection == 2) {
+   fKeepDfromB = kTRUE;
+   fKeepDfromBOnly = kFALSE;
+ }
+ return;
+
+}
+
+//______________________________________________________
+void AliCFVertexingHF::SetMCCandidateParam(Int_t label){
+       
+  fmcPartCandidate = dynamic_cast <AliAODMCParticle*> (fmcArray->At(label));
+  fNDaughters = fmcPartCandidate->GetNDaughters();
+  return;
+}
+
+
+//____________________________________________________________
+Int_t AliCFVertexingHF::MCcquarkCounting(AliAODMCParticle* mcPart) const{
+       
+  Int_t cquarks = 0;
+  if (mcPart->GetPdgCode() == 4) cquarks++; 
+  if (mcPart->GetPdgCode() == -4) cquarks++; 
+  if (!mcPart) {
+    AliWarning("Particle not found in tree, skipping\n"); 
+    return cquarks;
+  } 
+  
+  return cquarks;
+}
+
+
+//________________________________________________________
+Bool_t AliCFVertexingHF::CheckMCPartFamily(AliAODMCParticle */*mcPart*/, TClonesArray */*mcArray*/) const {
+
+  Int_t pdgGranma = CheckOrigin();
+  if (pdgGranma == -9999){
+    printf ("This particle come from a B decay channel but the we keep only the prompt charm particles\n");    
+    return kFALSE;
+  }    
+  
+ if (pdgGranma == -999){
+    printf ("This particle come from a prompt charm particles but we want only the ones coming from B\n");     
+    return kFALSE;
+  }    
+
+  if (!CheckMCDaughters()) return kFALSE;
+  if (!CheckMCChannelDecay()) return kFALSE;
+  return kTRUE;
+}
+
+//_________________________________________________________________________________________________
+Int_t AliCFVertexingHF::CheckOrigin() const {          
+  //
+  // checking whether the mother of the particles come from a charm or a bottom quark
+  //
+  
+  Int_t pdgGranma = 0;
+  Int_t mother = 0;
+  mother = fmcPartCandidate->GetMother();
+  Int_t istep = 0;
+  while (mother >0 ){
+    istep++;
+    AliDebug(2,Form("mother at step %d = %d", istep, mother));
+    AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
+    pdgGranma = mcGranma->GetPdgCode();
+    AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
+    Int_t abspdgGranma = TMath::Abs(pdgGranma);
+    if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
+      if (!fKeepDfromB) return -9999; //skip particle if come from a B meson.
+      
+      else{
+       return pdgGranma;
+      }
+    }
+    else {
+    if (fKeepDfromBOnly) return -999;
+    }
+    
+    mother = mcGranma->GetMother();
+  }
+  
+  return pdgGranma;
+}
+
+
+//___________________________________________
+Bool_t AliCFVertexingHF::CheckMCDaughters()const {
+       
+  AliAODMCParticle *mcPartDaughter;
+  Bool_t checkDaughters = kFALSE;
+
+  Int_t label0 = fmcPartCandidate->GetDaughter(0);
+  Int_t label1 = fmcPartCandidate->GetDaughter(1);
+  if (label1==0 || label0 == 0){
+    AliDebug(2, Form("The MC particle doesn't jave correct daughters, skipping!!"));
+       return checkDaughters;  
+  }
+
+  if (label1 - label0 != fProngs-1){
+    AliDebug(2, Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
+    return checkDaughters;  
+  }
+
+  for (Int_t iProng = 0; iProng<fProngs; iProng++){
+    mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));    
+    if (!mcPartDaughter) {
+      AliWarning("At least one Daughter Particle not found in tree, skipping"); 
+      return checkDaughters;  
+    }
+  }
+
+  checkDaughters = kTRUE;
+  return checkDaughters;
+}
+
+       
+
+//______________________________________________________
+Bool_t AliCFVertexingHF::FillMCContainer(Double_t *containerInputMC)
+{
+// fill the container for Generator level selection
+  Bool_t mcContainerFilled = kFALSE;  
+
+  Double_t* vectorMC = new Double_t[fNVar];
+  for (Int_t iVar = 0; iVar<fNVar; iVar++) vectorMC[iVar]= 9999.;
+  
+  if (GetGeneratedValuesFromMCParticle(&vectorMC[0])){
+    for (Int_t iVar = 0; iVar<fNVar; iVar++){
+      
+      containerInputMC[iVar] = vectorMC[iVar];
+    }
+    
+    mcContainerFilled = kTRUE;         
+  }
+  delete vectorMC;
+  return mcContainerFilled;    
+}
+
+//______________________________________________________
+Bool_t AliCFVertexingHF::FillRecoContainer(Double_t *containerInput) {  
+
+// fill the container for Reconstrucred level selection
+
+  Bool_t recoContainerFilled = kFALSE;
+  Double_t* vectorValues = new Double_t[fNVar];
+  Double_t* vectorReco = new Double_t[fNVar];  
+  
+  for (Int_t iVar = 0; iVar<fNVar; iVar++) {
+    vectorValues[iVar]= 9999.;
+    vectorReco[iVar]=9999.;
+  }
+
+  if(fFillFromGenerated){
+    //filled with MC values
+    if (GetGeneratedValuesFromMCParticle(&vectorValues[0])){
+      for (Int_t iVar = 0; iVar<fNVar; iVar++){
+                 containerInput[iVar] = vectorValues[iVar];
+               }
+               recoContainerFilled = kTRUE;            
+               }
+       }
+       else{
+    //filled with Reco values
+               
+               if (GetRecoValuesFromCandidate(&vectorReco[0])){
+                       for (Int_t iVar = 0; iVar<fNVar; iVar++){
+                               containerInput[iVar] = vectorReco[iVar];
+                       }
+                       recoContainerFilled = kTRUE;            
+               }
+       }
+  
+  delete vectorValues;
+  delete vectorReco;
+  return recoContainerFilled;  
+}
+
+//_____________________________________________________
+Bool_t AliCFVertexingHF::MCAcceptanceStep() const
+{
+
+  Bool_t bMCAccStep = kFALSE;
+  
+  AliAODMCParticle *mcPartDaughter;
+  Int_t label0 = fmcPartCandidate->GetDaughter(0);
+  Int_t label1 = fmcPartCandidate->GetDaughter(1);
+  if (label1==0 || label0 == 0){
+    AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
+    return bMCAccStep;  
+  }
+
+  if (label1 - label0 != fProngs-1){
+    AliDebug(2, Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
+       return bMCAccStep;  
+  }
+
+  for (Int_t iProng = 0; iProng<fProngs; iProng++){
+    mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));    
+    if (!mcPartDaughter) {
+      AliWarning("At least one Daughter Particle not found in tree, skipping"); 
+      return bMCAccStep;  
+    }
+    Double_t eta = mcPartDaughter->Eta();
+    Double_t pt = mcPartDaughter->Pt();
+         
+         //set values of eta and pt in the constructor.
+    if (TMath::Abs(eta) > 0.9 || pt < 0.1){
+               AliDebug(3,"At least one daughter has eta>0.9 or pt < 0.1 \n"); 
+      return bMCAccStep;
+    }
+
+  }
+
+  
+  bMCAccStep = kTRUE;
+  return bMCAccStep; 
+  
+}
+
+//_____________________________________________________
+Bool_t AliCFVertexingHF::MCRefitStep(AliAODEvent *aodEvent, AliESDtrackCuts *trackCuts) const
+{              
+  // check on the kTPCrefit and kITSrefit conditions of the daughters
+       Bool_t bRefitStep = kFALSE;
+
+       Int_t label0 = fmcPartCandidate->GetDaughter(0);
+       Int_t label1 = fmcPartCandidate->GetDaughter(1);
+
+       if (label1==0 || label0 == 0){
+               AliDebug(2, Form("The MC particle doesn't jave correct daughters, skipping!!"));
+        return bRefitStep;  
+       }
+
+       if (label1 - label0 != fProngs-1){
+               AliDebug(2, Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
+               //AliInfo(Form("The MC particle doesn't come from a %d-prong decay, skipping!!", fProngs));
+               return bRefitStep;  
+       }
+
+       Int_t foundDaughters = 0;
+
+       if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
+    
+               for(Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
+                       AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
+                       if (track->GetLabel()>= label0 && track->GetLabel()<=label1){
+                               foundDaughters++;
+                               printf("daughter %d \n",foundDaughters);
+                               if(trackCuts->GetRequireTPCRefit()){
+                                       if((track->GetStatus()&AliESDtrack::kTPCrefit)) {
+                                               bRefitStep = kTRUE;
+                                       }
+                                       else {
+                                               AliDebug(3, "Refit cut not passed , missing TPC refit\n");
+                                               AliInfo( "Refit cut not passed , missing TPC refit\n");
+                                               return kFALSE;
+                                       }
+                               }
+                               
+                               if (trackCuts->GetRequireITSRefit()) {
+                                       if((track->GetStatus()&AliESDtrack::kITSrefit)){
+                                               bRefitStep = kTRUE;
+                                       }
+                                       else {
+                                               AliDebug(3, "Refit cut not passed , missing ITS refit\n");
+                                               AliInfo("Refit cut not passed , missing ITS refit\n");
+                                               return kFALSE;
+                                       }
+                               }
+                       }       
+                       if (foundDaughters == fProngs){
+        
+                               break;
+                       }
+      
+               }
+    
+       }
+  
+       return bRefitStep;
+}
+
+//____________________________________________________________________________
+
+Bool_t AliCFVertexingHF::RecoStep() 
+//check also vertex and ITS Refit and TPC Refit
+{ 
+  Bool_t bRecoStep = kFALSE;
+  Int_t mcLabel = GetMCLabel();
+  
+  if (mcLabel == -1) {
+    AliDebug(2,"No MC particle found");
+    return bRecoStep;
+  }
+  else{
+    fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
+    if (!fmcPartCandidate){
+      AliWarning("Could not find associated MC in AOD MC tree");
+      return bRecoStep;
+    }
+    
+  }
+  
+  Int_t pdgGranma = CheckOrigin();
+  
+  if (pdgGranma == -9999){
+    printf ("This particle come from a B decay channel but we keep only prompt charm particles\n");
+    return bRecoStep;
+  }
+
+  if (pdgGranma == -999){
+    printf ("This particle come from a  prompt charm particle but we want only the ones coming from B\n");
+    return bRecoStep;
+  }
+  
+   
+  bRecoStep=kTRUE;
+  return bRecoStep;
+  
+  
+}      
+//____________________________________________
+Double_t AliCFVertexingHF::GetEtaProng(Int_t iProng) const {
+
+  if (fRecoCandidate){
+    Double_t etaProng = fRecoCandidate->EtaProng(iProng);  
+    return etaProng;
+  }
+  return 999999;  
+}
+//______________________________________________________
+Double_t AliCFVertexingHF::GetPtProng(Int_t iProng) const {
+
+  if (fRecoCandidate){
+    Double_t ptProng = fRecoCandidate->PtProng(iProng);  
+    return ptProng;
+  }
+  return 999999;  
+  
+}
+
+//____________________________________________________________________
+
+Bool_t AliCFVertexingHF::RecoAcceptStep(AliESDtrackCuts *trackCuts) const
+{
+  Bool_t bRecoAccStep = kFALSE;
+  
+  Float_t etaCutMin, ptCutMin, etaCutMax, ptCutMax;
+  trackCuts->GetEtaRange(etaCutMin, etaCutMax);
+  trackCuts->GetPtRange(ptCutMin, ptCutMax);
+  
+  Float_t etaProng=0., ptProng=0.; 
+  
+  for (Int_t iProng =0; iProng<fProngs; iProng++){
+    
+    etaProng = GetEtaProng(iProng);
+    ptProng = GetPtProng(iProng);
+    
+    Bool_t acceptanceProng = (etaProng>etaCutMin && etaProng<etaCutMax && ptProng>ptCutMin && ptProng<ptCutMax);
+    if (!acceptanceProng) {
+      printf ("At least one reconstructed prong isn't in the acceptance\n");
+      return bRecoAccStep;
+    }
+  }
+  
+  bRecoAccStep=kTRUE;
+  return bRecoAccStep;
+}
+//___________________________________________________________
+
+Bool_t AliCFVertexingHF::FillUnfoldingMatrix(Double_t *fill) const{
+
+  fill = new Double_t[4];
+  
+  if(fmcPartCandidate){
+    
+    fill[0] = GetPtCand();
+    fill[1] = GetYCand();
+    
+    fill[2] =  fmcPartCandidate->Pt(); 
+    fill[3] =  fmcPartCandidate->Y(); 
+    
+    return kTRUE;
+  }
+
+  delete fill;
+  return kFALSE;
+}
+
+
+//___________________________________________________________
diff --git a/PWG3/vertexingHF/AliCFVertexingHF.h b/PWG3/vertexingHF/AliCFVertexingHF.h
new file mode 100644 (file)
index 0000000..33eeda8
--- /dev/null
@@ -0,0 +1,115 @@
+#ifndef ALICFVERTEXINGHF_H
+#define ALICFVERTEXINGHF_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 and step 
+// Author : C. Zampolli, CERN
+// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
+// Base class for HF Unfolding - agrelli@uu.nl
+//-----------------------------------------------------------------------
+
+#include "AliCFContainer.h"
+#include "AliAODRecoDecayHF.h"
+
+class TH1I;
+class TParticle ;
+class TFile ;
+class TClonesArray ;
+class AliAODMCParticle;
+class AliAODMCHeader;
+class AliAODEvent;
+class THnSparse;
+class TClonesArray;
+class AliESDtrackCuts;
+
+
+
+class AliCFVertexingHF : public TObject {
+       public:
+       
+       AliCFVertexingHF() ;
+       AliCFVertexingHF(TClonesArray *mcArray, UShort_t originDselection);
+       AliCFVertexingHF(const AliCFVertexingHF& c);
+       AliCFVertexingHF& operator= (const AliCFVertexingHF& c);
+
+       virtual ~AliCFVertexingHF();
+               
+       virtual Bool_t GetGeneratedValuesFromMCParticle(Double_t* /*vectorMC*/) {return kFALSE;} 
+       virtual Bool_t GetRecoValuesFromCandidate(Double_t* /*vectorReco*/) const {return kFALSE;}
+       virtual Bool_t CheckMCChannelDecay() const {return kFALSE;}
+       
+       void   SetFillFromGenerated(Bool_t flag) {fFillFromGenerated = flag;}
+       Bool_t GetFillFromGenerated() const {return fFillFromGenerated;}
+               
+       void  SetNVar(Int_t nVar) {fNVar = nVar;}  
+
+       void  SetRecoPrimVertex (Double_t zPrimVertex) {fzPrimVertex = zPrimVertex;}
+       void  SetMCPrimaryVertex (Double_t zMCVertex){fzMCVertex = zMCVertex;}
+       void  SetMCLabel (Int_t mcLabel) {fmcLabel = mcLabel;}
+       Int_t GetMCLabel () const {return  fmcLabel;}
+               
+       void   SetMCCandidateParam(Int_t label);
+
+
+       Int_t  MCcquarkCounting(AliAODMCParticle* mcPart) const; 
+       Bool_t CheckMCPartFamily(AliAODMCParticle */*mcPart*/, TClonesArray */*mcArray*/) const;
+       //      Int_t  CheckOrigin(AliAODMCParticle* mcPart) const;
+       Int_t  CheckOrigin() const;
+       Bool_t CheckMCDaughters() const;
+       Bool_t FillMCContainer(Double_t *containerInputMC); 
+       Bool_t FillRecoContainer(Double_t *containerInput); 
+       Bool_t MCAcceptanceStep() const;
+       Bool_t MCRefitStep(AliAODEvent *aodEvent, AliESDtrackCuts *trackCuts) const;
+       Bool_t RecoStep();
+
+       Double_t GetEtaProng(Int_t iProng) const;
+       Double_t GetPtProng(Int_t iProng) const;
+
+       Double_t GetPtCand() const {return fRecoCandidate->Pt();}
+       Double_t GetYCand() const {return fRecoCandidate->Y();}
+
+       Bool_t RecoAcceptStep(AliESDtrackCuts *trackCuts) const;
+       
+       Bool_t FillUnfoldingMatrix(Double_t *fill) const;
+       
+       void SetNProngs(Int_t nProngs){fProngs = nProngs;}
+       void SetDselection(UShort_t originDselection); 
+
+       protected:
+       
+       TClonesArray      *fmcArray;               //mcArray candidate
+       AliAODRecoDecayHF *fRecoCandidate;         // Reconstructed HF candidate 
+       AliAODMCParticle  *fmcPartCandidate;
+       Int_t fNDaughters;
+       Int_t fNVar;                // get Number of variables for the container from the channel decay
+       Double_t fzPrimVertex;       //Reco z primary vertex    
+       Double_t fzMCVertex;         //MC z primary vertex
+       
+       Bool_t fFillFromGenerated;   //  flag to indicate whether data container should be filled  
+       UShort_t fOriginDselection;      // flag to select D0 origins. 0 Only from charm 1 only from beauty 2 both from charm and beauty
+       
+       Bool_t fKeepDfromB;           //flag for the feed down from b quark decay.                      
+       Bool_t fKeepDfromBOnly;       // flag to keep only the charm particles that comes from beauty decays
+       Int_t fmcLabel;              // results of the MatchToMC()
+       Int_t fProngs;               // n. of prongs    
+       ClassDef(AliCFVertexingHF, 1);
+       
+};
+
+#endif
diff --git a/PWG3/vertexingHF/AliCFVertexingHF2Prong.cxx b/PWG3/vertexingHF/AliCFVertexingHF2Prong.cxx
new file mode 100644 (file)
index 0000000..596554e
--- /dev/null
@@ -0,0 +1,292 @@
+/**************************************************************************
+ * 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 and step 
+// Author : C. Zampolli, CERN
+// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
+// Base class for HF Unfolding - agrelli@uu.nl
+//-----------------------------------------------------------------------
+
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODMCParticle.h"
+#include "AliAODEvent.h"
+#include "TClonesArray.h"
+#include "AliCFVertexingHF.h"
+#include "AliESDtrack.h"
+#include "TDatabasePDG.h"
+
+#include "AliCFVertexingHF2Prong.h"
+#include "AliCFContainer.h"
+
+ClassImp(AliCFVertexingHF2Prong)
+
+
+//_________________________________________
+  AliCFVertexingHF2Prong::AliCFVertexingHF2Prong(TClonesArray *mcArray, UShort_t originDselection):
+  AliCFVertexingHF(mcArray, originDselection)
+{      
+       
+  SetNProngs(2);
+}
+
+
+//_____________________________________
+AliCFVertexingHF2Prong& AliCFVertexingHF2Prong::operator=(const AliCFVertexingHF2Prong& c){
+  if  (this != &c) {
+
+         AliCFVertexingHF::operator=(c);
+   
+  }
+    return *this;
+
+
+}
+
+//__________________________________________
+Bool_t AliCFVertexingHF2Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand){
+  
+  Bool_t bSignAssoc = kFALSE;
+
+  fRecoCandidate = recoCand;
+
+  if (!fRecoCandidate) {
+    printf("fRecoCandidate not found, problem in assignement\n");
+    return bSignAssoc;
+  }
+  
+  if (fRecoCandidate->GetPrimaryVtx()) printf("fReco Candidate has a pointer to PrimVtx\n");
+  if (recoCand->GetPrimaryVtx()) printf("Reco Cand has a pointer to PrimVtx\n");
+  
+  Int_t pdgCand = 421;
+  Int_t pdgDgD0toKpi[2]={321,211};
+  Int_t nentries = fmcArray->GetEntriesFast();
+
+  printf ("nentries = %d\n", nentries);
+  Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,2,pdgDgD0toKpi);
+  
+  if (mcLabel == -1) return bSignAssoc;
+  SetMCLabel(mcLabel);
+  fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
+    
+  if (!fmcPartCandidate){
+    AliDebug(3,"No part candidate");
+    return bSignAssoc;
+  }
+
+
+  bSignAssoc = kTRUE;
+  return bSignAssoc;
+}
+
+//______________________________________________
+Bool_t AliCFVertexingHF2Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) {
+       // 
+       // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
+       //
+       
+       Bool_t bGenValues = kFALSE;
+       
+       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
+       fmcPartCandidate->XvYvZv(vtx1);  // cm
+
+       
+       Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
+       Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
+       AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
+       AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
+
+       // 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 bGenValues;
+       }
+       
+       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 (fmcPartCandidate->GetPdgCode() == 421){  // D0
+         AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
+         cosThetaStar = cosThetaStarD0;
+       }
+       else if (fmcPartCandidate->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 bGenValues;
+       }
+               
+       //ct
+       Double_t cT = decay->Ct(421);
+       // 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] = fmcPartCandidate->Pt();
+       vectorMC[1] = fmcPartCandidate->Y() ;
+       vectorMC[2] = cosThetaStar ;
+       vectorMC[3] = pTpi ;
+       vectorMC[4] = pTK ;
+       vectorMC[5] = cT*1.E4 ;  // in micron
+       vectorMC[6] = 0.;   // dummy value, meaningless in MC
+       vectorMC[7] = 0.;   // dummy value, meaningless in MC, in micron
+       vectorMC[8] = 0.;   // dummy value, meaningless in MC, in micron
+       vectorMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
+       vectorMC[10] = 1.01;    // dummy value, meaningless in MC
+       vectorMC[11] = fmcPartCandidate->Phi(); 
+       vectorMC[12] = fzMCVertex;    // z of reconstructed of primary vertex
+       bGenValues = kTRUE;
+       return bGenValues;
+}
+
+
+
+//____________________________________________
+Bool_t AliCFVertexingHF2Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
+{ 
+  Bool_t bFillRecoValues=kFALSE;
+  
+  AliAODRecoDecayHF2Prong *d0toKpi = (AliAODRecoDecayHF2Prong*)fRecoCandidate;
+  if (d0toKpi->GetPrimaryVtx())printf("d0toKpi has primary vtx\n");
+  if (fRecoCandidate->GetPrimaryVtx())printf("fRecoCandidate has primary vtx\n");
+  
+  Double_t pt = d0toKpi->Pt();
+  Double_t rapidity = d0toKpi->YD0();
+  Double_t invMass=0.;
+  Double_t cosThetaStar = 9999.;
+  Double_t pTpi = 0.;
+  Double_t pTK = 0.;
+  Double_t dca = d0toKpi->GetDCA();
+  Double_t d0pi = 0.;
+  Double_t d0K = 0.;
+  Double_t d0xd0 = d0toKpi->Prodd0d0();
+  Double_t cosPointingAngle = d0toKpi->CosPointingAngle();
+  Double_t phi = d0toKpi->Phi();
+  Int_t pdgCode = fmcPartCandidate->GetPdgCode();
+
+  if (pdgCode > 0){
+    cosThetaStar = d0toKpi->CosThetaStarD0();
+    pTpi = d0toKpi->PtProng(0);
+    pTK = d0toKpi->PtProng(1);
+    d0pi = d0toKpi->Getd0Prong(0);
+    d0K = d0toKpi->Getd0Prong(1);
+    invMass=d0toKpi->InvMassD0();
+  }
+  else {
+    cosThetaStar = d0toKpi->CosThetaStarD0bar();
+    pTpi = d0toKpi->PtProng(1);
+    pTK = d0toKpi->PtProng(0);
+    d0pi = d0toKpi->Getd0Prong(1);
+    d0K = d0toKpi->Getd0Prong(0);
+    invMass= d0toKpi->InvMassD0bar();
+  }
+  
+  Double_t cT = d0toKpi->CtD0();
+  
+  
+  vectorReco[0] = pt;
+  vectorReco[1] = rapidity;
+  vectorReco[2] = cosThetaStar;
+  vectorReco[3] = pTpi;
+  vectorReco[4] = pTK;
+  vectorReco[5] = cT*1.E4;  // in micron
+  vectorReco[6] = dca*1.E4;  // in micron
+  vectorReco[7] = d0pi*1.E4;  // in micron
+  vectorReco[8] = d0K*1.E4;  // in micron
+  vectorReco[9] = d0xd0*1.E8;  // in micron^2
+  vectorReco[10] = cosPointingAngle;  // in micron
+  vectorReco[11] = phi;  
+  vectorReco[12] = fzPrimVertex;    // z of reconstructed of primary vertex
+  
+  bFillRecoValues = kTRUE;
+
+
+  return bFillRecoValues;
+}
+
+
+//_____________________________________________________________
+Bool_t AliCFVertexingHF2Prong::CheckMCChannelDecay() const
+{ 
+  Bool_t checkCD = kFALSE;
+  Int_t daughter0 = fmcPartCandidate->GetDaughter(0);
+  Int_t daughter1 = fmcPartCandidate->GetDaughter(1);
+  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter0));
+  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter1));
+
+  if (!mcPartDaughter0 || !mcPartDaughter1) {
+    AliDebug (2,"Problems in the MC Daughters\n");
+    return checkCD;
+  }
+
+  if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
+       TMath::Abs(mcPartDaughter1->GetPdgCode())==211) && 
+      !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
+       TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
+    AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
+    return checkCD;  
+  }
+  
+  checkCD = kTRUE;
+  return checkCD;
+  
+}
+
+//_______________________________________________
diff --git a/PWG3/vertexingHF/AliCFVertexingHF2Prong.h b/PWG3/vertexingHF/AliCFVertexingHF2Prong.h
new file mode 100644 (file)
index 0000000..f69c9ae
--- /dev/null
@@ -0,0 +1,63 @@
+#ifndef ALICFVERTEXINGHF2PRONG_H
+#define ALICFVERTEXINGHF2PRONG_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 and step 
+// Author : C. Zampolli, CERN
+// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
+// Base class for HF Unfolding - agrelli@uu.nl
+//-----------------------------------------------------------------------
+
+
+#include "AliCFVertexingHF.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoDecayHF2Prong.h"
+
+class AliAODMCParticle;
+class TClonesArray;
+class AliCFVertexingHF;
+class AliESDtrack;
+class TDatabasePDG;
+
+class AliCFVertexingHF2Prong : public AliCFVertexingHF{
+       public:
+               
+  AliCFVertexingHF2Prong(){};
+  AliCFVertexingHF2Prong(TClonesArray *mcArray, UShort_t originDselection);
+       
+  virtual ~AliCFVertexingHF2Prong(){};
+  
+  
+  Bool_t GetGeneratedValuesFromMCParticle(Double_t* /*vectorMC*/);
+  Bool_t GetRecoValuesFromCandidate(Double_t* /*vectorReco*/ ) const;
+  Bool_t CheckMCChannelDecay()const;
+  
+  Bool_t SetRecoCandidateParam(AliAODRecoDecayHF *recoCand);
+  
+ protected:
+  
+  
+ private:      
+  AliCFVertexingHF2Prong(const AliCFVertexingHF2Prong& c);
+  AliCFVertexingHF2Prong& operator= (const AliCFVertexingHF2Prong& other);
+  
+  ClassDef(AliCFVertexingHF2Prong, 1);
+  
+};
+
+#endif
diff --git a/PWG3/vertexingHF/macros/AddTaskCFVertexingHF.C b/PWG3/vertexingHF/macros/AddTaskCFVertexingHF.C
new file mode 100644 (file)
index 0000000..008077a
--- /dev/null
@@ -0,0 +1,500 @@
+//DEFINITION OF A FEW CONSTANTS
+const Double_t ymin  = -2.1 ;
+const Double_t ymax  =  2.1 ;
+// 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.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 Double_t phimax = 2Pi;  // defined in the macro!!!!!!!!!!!!!!  
+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 ;
+// cuts
+const Double_t ptmin = 0.1;
+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 = 5;
+
+//----------------------------------------------------
+
+AliCFTaskVertexingHF *AddTaskCFVertexingHF(const char* cutFile = "./D0toKpiCuts.root",Bool_t isKeepDfromB=kFALSE, Bool_t isKeepDfromBOnly=kFALSE)
+{
+       printf("Addig CF task using cuts from file %s\n",cutFile);
+
+       TFile* fileCuts = new TFile(cutFile);
+       AliRDHFCutsD0toKpi *cutsD0toKpi = (AliRDHFCutsD0toKpi*)fileCuts->Get("D0toKpiCuts");
+       
+       // check that the fKeepD0fromB flag is set to true when the fKeepD0fromBOnly flag is true
+       //  for now the binning is the same than for all D's
+       if(isKeepDfromBOnly) isKeepDfromB = true;
+       
+
+       /*
+         Double_t ptmin_0_4;
+         Double_t ptmax_0_4;
+         Double_t ptmin_4_8;
+         Double_t ptmax_4_8;
+         Double_t ptmin_8_10;
+         Double_t ptmax_8_10;
+         
+         if(!isKeepDfromB){
+         ptmin_0_4 =  0.0 ;
+         ptmax_0_4 =  4.0 ;
+         ptmin_4_8 =  4.0 ;
+         ptmax_4_8 =  8.0 ;
+         ptmin_8_10 =  8.0 ;
+         ptmax_8_10 =  10.0 ;
+         } else{
+         ptmin_0_4 =  0.0 ;
+         ptmax_0_4 =  3.0 ;
+         ptmin_4_8 =  3.0 ;
+         ptmax_4_8 =  5.0 ;
+         ptmin_8_10 =  5.0 ;
+         ptmax_8_10 =  10.0 ;
+         }
+       */
+
+       //CONTAINER DEFINITION
+       Info("AliCFTaskVertexingHF","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 ipTk  = 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;
+
+       const Double_t phimax = 2*TMath::Pi();
+
+       //Setting up the container grid... 
+       UInt_t nstep = 10; //number of selection steps: MC with limited acceptance, MC, Acceptance, Vertex, Refit, Reco (no cuts), RecoAcceptance, RecoITSClusters (RecoAcceptance included), RecoPPR (RecoAcceptance+RecoITSCluster included), RecoPID 
+       const Int_t nvar   = 13 ; //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  = 1 ; //bins in pt from 8 to 10 GeV
+
+/*
+       Int_t nbin0_0_4;
+       Int_t nbin0_4_8;
+       Int_t nbin0_8_10;
+       if (!isKeepDfromB){
+         nbin0_0_4  = 8 ; //bins in pt from 0 to 4 GeV
+         nbin0_4_8  = 4 ; //bins in pt from 4 to 8 GeV
+         nbin0_8_10  = 1 ; //bins in pt from 8 to 10 GeV
+       }else{
+         nbin0_0_4  = 3 ; //bins in pt from 0 to 3 GeV
+         nbin0_4_8  = 1 ; //bins in pt from 3 to 5 GeV
+         nbin0_8_10  = 1 ; //bins in pt from 5 to 10 GeV
+       }
+*/
+       const Int_t nbin0 = cutsD0toKpi->GetNPtBins(); // bins in pT
+       printf("pT: nbin (from cuts file) = %d\n",nbin0);
+       const Int_t nbin1  = 42 ; //bins in y
+       const Int_t nbin2  = 42 ; //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
+       const Int_t nbin6  = 24 ; //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  = 1050 ; //bins in cosPointingAngle
+       const Int_t nbin11  = 20 ; //bins in Phi
+       const Int_t nbin12  = 60 ; //bins in z vertex
+
+       //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[0]=nbin0;
+       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[3]=nbin0;
+       iBin[4]=nbin0;
+       iBin[5]=nbin5;
+       iBin[6]=nbin6;
+       iBin[7]=nbin7;
+       iBin[8]=nbin8;
+       iBin[9]=nbin9;
+       iBin[10]=nbin10;
+       iBin[11]=nbin11;
+       iBin[12]=nbin12;
+       
+       //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];
+
+       // checking limits
+       /*
+       if (ptmax_0_4 != ptmin_4_8) {
+               Error("AliCFHeavyFlavourTaskMultiVarMultiStep","max lim 1st range != min lim 2nd range, please check!");
+       }
+       if (ptmax_4_8 != ptmin_8_10) {
+               Error("AliCFHeavyFlavourTaskMultiVarMultiStep","max lim 2nd range != min lim 3rd range, please check!");
+       }
+       */
+       // values for bin lower bounds
+       // pt
+       Float_t* floatbinLim0 = cutsD0toKpi->GetPtBinLimits();
+       for (Int_t ibin0 = 0 ; ibin0<iBin[0]+1; ibin0++){
+               binLim0[ibin0] = (Double_t)floatbinLim0[ibin0];
+               binLim3[ibin0] = (Double_t)floatbinLim0[ibin0];
+               binLim4[ibin0] = (Double_t)floatbinLim0[ibin0];
+       }
+       for(Int_t i=0; i<=nbin0; i++) printf("binLim0[%d]=%f\n",i,binLim0[i]);  
+
+       /*
+       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("AliCFHeavyFlavourTaskMultiVarMultiStep","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("AliCFHeavyFlavourTaskMultiVarMultiStep","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("AliCFHeavyFlavourTaskMultiVarMultiStep","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("AliCFHeavyFlavourTaskMultiVarMultiStep","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("AliCFHeavyFlavourTaskMultiVarMultiStep","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("AliCFHeavyFlavourTaskMultiVarMultiStep","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 ;
+
+       // 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 ;
+               //              Info("AliCFHeavyFlavourTaskMultiVarMultiStep",Form("i-th bin, lower limit = %f", binLim12[i]));
+       }
+
+       // debugging printings
+       //Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Printing lower limits for bins in pt");
+       //for (Int_t i =0; i<= iBin[0]; i++){
+       //      Info("AliCFHeavyFlavourTaskMultiVarMultiStep",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("AliCFHeavyFlavourTaskMultiVarMultiStep",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("AliCFHeavyFlavourTaskMultiVarMultiStep",Form("i-th bin, lower limit = %f", binLim4[i]));
+       //      }
+
+       //one "container" for MC
+       TString nameContainer="";
+       if(!isKeepDfromB) {
+         nameContainer="CFHFccontainer0_New";
+       }
+       else  if(isKeepDfromBOnly){
+         nameContainer="CFHFccontainer0DfromB_New";
+       }
+       else  {
+         nameContainer="CFHFccontainer0allD_New";
+         
+       }
+       
+       AliCFContainer* container = new AliCFContainer(nameContainer,"container for tracks",nstep,nvar,iBin);
+       //setting the bin limits
+       printf("pt\n");
+       container -> SetBinLimits(ipt,binLim0);
+       printf("y\n");
+       container -> SetBinLimits(iy,binLim1);
+       printf("cts\n");
+       container -> SetBinLimits(icosThetaStar,binLim2);
+       printf("ptPi\n");
+       container -> SetBinLimits(ipTpi,binLim3);
+       printf("ptK\n");
+       container -> SetBinLimits(ipTk,binLim4);
+       printf("cT\n");
+       container -> SetBinLimits(icT,binLim5);
+       printf("dca\n");
+       container -> SetBinLimits(idca,binLim6);
+       printf("d0Pi\n");
+       container -> SetBinLimits(id0pi,binLim7);
+       printf("d0K\n");
+       container -> SetBinLimits(id0K,binLim8);
+       printf("d0xd0\n");
+       container -> SetBinLimits(id0xd0,binLim9);
+       printf("pointing\n");
+       container -> SetBinLimits(ipointing,binLim10);
+       printf("phi\n");
+       container -> SetBinLimits(iphi,binLim11);
+       printf("z\n");
+       container -> SetBinLimits(iz,binLim12);
+       
+       container -> SetStepTitle(0, "MCLimAcc");
+       container -> SetStepTitle(1, "MC");
+        container -> SetStepTitle(2, "MCAcc");
+        container -> SetStepTitle(3, "RecoVertex");
+        container -> SetStepTitle(4, "RecoRefit");
+        container -> SetStepTitle(5, "Reco");
+        container -> SetStepTitle(6, "RecoAcc");
+       container -> SetStepTitle(7, "RecoITSCluster");
+       container -> SetStepTitle(8, "RecoCuts");
+       container -> SetStepTitle(8, "RecoPID");
+
+        container -> SetVarTitle(ipt,"pt");
+       container -> SetVarTitle(iy,"y");
+        container -> SetVarTitle(icosThetaStar, "cosThetaStar");
+        container -> SetVarTitle(ipTpi, "ptpi");
+       container -> SetVarTitle(ipTk, "ptK");
+        container -> SetVarTitle(icT, "ct");
+        container -> SetVarTitle(idca, "dca");
+        container -> SetVarTitle(id0pi, "d0pi");
+        container -> SetVarTitle(id0K, "d0K");
+       container -> SetVarTitle(id0xd0, "d0xd0");
+       container -> SetVarTitle(ipointing, "piointing");
+       container -> SetVarTitle(iphi, "phi");
+       container -> SetVarTitle(iz, "z");
+
+
+       //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(PDG, kTRUE);  // kTRUE set in order to include D0_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, Limited Acceptance
+       man->SetParticleCutsList(1 , mcList); // MC
+       man->SetParticleCutsList(2 , accList); // Acceptance 
+       man->SetParticleCutsList(3 , emptyList); // Vertex 
+       man->SetParticleCutsList(4 , emptyList); // Refit 
+       man->SetParticleCutsList(5 , emptyList); // AOD
+       man->SetParticleCutsList(6 , emptyList); // AOD in Acceptance
+       man->SetParticleCutsList(7 , emptyList); // AOD with required n. of ITS clusters
+       man->SetParticleCutsList(8 , emptyList); // AOD Reco (PPR cuts implemented in Task)
+       man->SetParticleCutsList(9 , emptyList); // AOD Reco PID
+       
+       // 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
+       AliCFTaskVertexingHF *task = new AliCFTaskVertexingHF("AliCFTaskVertexingHF",cutsD0toKpi);
+       // task->SetFillFromGenerated(kFALSE);
+       //task->SetMinITSClusters(minITSClusters);
+       task->SetCFManager(man); //here is set the CF manager
+       //      task->SetFeedDownExclusion(isKeepDfromB);
+       
+
+        //-----------------------------------------------------------//
+        //   create correlation matrix for unfolding - only eta-pt   //
+        //-----------------------------------------------------------//
+
+        Bool_t AcceptanceUnf = kTRUE; // unfold at acceptance level, otherwise PPR
+
+        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];
+
+       TString nameCorr="";
+       if(!isKeepDfromB) {
+         nameCorr="CFHFcorr0_New";
+       }
+       else  if(isKeepD0fromBOnly){
+         nameCorr= "CFHFcorr0KeepDfromBOnly";
+       }
+       else  {
+         nameCorr="CFHFcorr0allD_New";
+
+       }
+
+        THnSparseD* correlation = new THnSparseD(nameCorr,"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();
+       TString output1name="", output2name="", output3name="";
+       output2name=nameContainer;
+       output3name=nameCorr;
+       if(!isKeepDfromB) {
+         outputfile += ":PWG3_D2H_CFtaskD0toKpi_NEW";
+         output1name="CFHFchist0_New";
+       }
+       else  if(isKeepDfromBOnly){
+         outputfile += ":PWG3_D2H_CFtaskD0toKpiKeepDfromBOnly";
+         output1name="CFHFchist0DfromB";
+       }
+       else{
+         outputfile += ":PWG3_D2H_CFtaskD0toKpiKeepDfromB_NEW";
+         output1name="CFHFchist0allD_New";
+         
+       }
+
+       //now comes user's output objects :
+       // output TH1I for event counting
+       AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(output1name, TH1I::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+       // output Correction Framework Container (for acceptance & efficiency calculations)
+       AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(output2name, AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+       // Unfolding - correlation matrix
+        AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(output3name, 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;
+}
+