New class to handle multi-dimensional unfolding + user macros (test/)
[u/mrichter/AliRoot.git] / CORRFW / test / AliCFTaskForUnfolding.C
1 //DEFINITION OF A FEW CONSTANTS
2 const Double_t ptmin  =   0.0 ;
3 const Double_t ptmax  =   1.0 ;
4 const Double_t etamin =  -1.5 ;
5 const Double_t etamax =   1.5 ;
6 const Int_t    PDG    =    211; 
7 const Int_t    minclustersTPC = 40 ;
8 //----------------------------------------------------
9
10 Bool_t AliCFTaskForUnfolding()
11 {
12   
13   TBenchmark benchmark;
14   benchmark.Start("AliSingleTrackTask");
15
16   AliLog::SetGlobalDebugLevel(0);
17
18   Load() ; //load the required libraries
19
20   TChain * analysisChain ;
21   printf("\n\nRunning on local file, please check the path\n\n");
22
23   analysisChain = new TChain("esdTree");
24   analysisChain->Add("AliESDs.root");
25   
26   Info("AliCFTaskForUnfolding",Form("CHAIN HAS %d ENTRIES",(Int_t)analysisChain->GetEntries()));
27
28   //CONTAINER DEFINITION
29   Info("AliCFTaskForUnfolding","SETUP CONTAINER");
30   //the sensitive variables (2 in this example), their indices
31   UInt_t ipt = 0;
32   UInt_t ieta  = 1;
33   //Setting up the container grid... 
34   UInt_t nstep = 3 ; //number of selection steps MC 
35   const Int_t nvar   = 2 ; //number of variables on the grid:pt,eta
36   const Int_t nbin[nvar] = {20,20} ;
37
38   //arrays for the number of bins in each dimension
39   Int_t iBin[nvar];
40   iBin[0]=nbin[0];
41   iBin[1]=nbin[1];
42
43   //arrays for lower bounds :
44   Double_t *binLim0=new Double_t[nbin[0]+1];
45   Double_t *binLim1=new Double_t[nbin[1]+1];
46
47   //values for bin lower bounds
48   for(Int_t i=0; i<=nbin[0]; i++) binLim0[i]=(Double_t)ptmin  + ( ptmax- ptmin)/nbin[0]*(Double_t)i ; 
49   for(Int_t i=0; i<=nbin[1]; i++) binLim1[i]=(Double_t)etamin + (etamax-etamin)/nbin[1]*(Double_t)i ; 
50
51   //one "container" for MC
52   AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin);
53   //setting the bin limits
54   container -> SetBinLimits(ipt,binLim0);
55   container -> SetBinLimits(ieta ,binLim1);
56
57   // Gen-Level kinematic cuts
58   AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
59   mcKineCuts->SetPtRange (ptmin ,ptmax );
60   mcKineCuts->SetEtaRange(etamin,etamax);
61
62   //Particle-Level cuts:  
63   AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
64   mcGenCuts->SetRequireIsPrimary();
65   mcGenCuts->SetRequirePdgCode(PDG);
66
67   // Rec-Level kinematic cuts
68   AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");
69   recKineCuts->SetPtRange( ptmin, ptmax);
70   recKineCuts->SetPtRange(etamin,etamax);
71
72   AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
73   recQualityCuts->SetMinNClusterTPC(minclustersTPC);
74
75   AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
76   recIsPrimaryCuts->SetMaxNSigmaToVertex(3);
77
78   printf("CREATE MC KINE CUTS\n");
79   TObjArray* mcList = new TObjArray(0) ;
80   mcList->AddLast(mcKineCuts);
81   mcList->AddLast(mcGenCuts);
82
83   printf("CREATE RECONSTRUCTION CUTS\n");
84   TObjArray* recList = new TObjArray(0) ;
85   recList->AddLast(recKineCuts);
86   recList->AddLast(recQualityCuts);
87   recList->AddLast(recIsPrimaryCuts);
88
89   TObjArray* emptyList = new TObjArray(0);
90
91   //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
92   printf("CREATE INTERFACE AND CUTS\n");
93   AliCFManager* man = new AliCFManager() ;
94   man->SetNStepEvent(0);
95   man->SetParticleContainer(container);
96   man->SetParticleCutsList(0,mcList);
97   man->SetParticleCutsList(1,recList);
98   man->SetParticleCutsList(2,emptyList); // this is special step for monte carlo spectrum
99
100   //CREATE THE TASK
101   printf("CREATE TASK\n");
102   // create the task
103   AliCFTaskForUnfolding *task = new AliCFTaskForUnfolding("AliCFTaskForUnfolding");
104   task->SetCFManager(man); //here is set the CF manager
105
106   //create correlation matrix for unfolding
107   Int_t thnDim[2*nvar];
108   for (int k=0; k<nvar; k++) {
109     //first half  : reconstructed 
110     //second half : MC
111     thnDim[k]      = nbin[k];
112     thnDim[k+nvar] = nbin[k];
113   }
114   THnSparseD* correlation = new THnSparseD("correlation","THnSparse with correlations",2*nvar,thnDim);
115   Double_t** binEdges = new Double_t[nvar];
116   for (int k=0; k<nvar; k++) {
117     binEdges[k]=new Double_t[nbin[k]+1];
118     container->GetBinLimits(k,binEdges[k]);
119     correlation->SetBinEdges(k,binEdges[k]);
120     correlation->SetBinEdges(k+nvar,binEdges[k]);
121   }
122   task->SetCorrelationMatrix(correlation);
123   //---
124
125   //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS
126   printf("CREATE ANALYSIS MANAGER\n");
127   // Make the analysis manager
128   AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
129   mgr->SetAnalysisType(AliAnalysisManager::kLocalAnalysis);
130
131   AliMCEventHandler*  mcHandler = new AliMCEventHandler();
132   mgr->SetMCtruthEventHandler(mcHandler);
133  
134   AliInputEventHandler* dataHandler = new AliESDInputHandler();
135   mgr->SetInputEventHandler(dataHandler);
136
137   // Create and connect containers for input/output
138
139   //------ input data ------
140   AliAnalysisDataContainer *cinput0  = mgr->CreateContainer("cchain0",TChain::Class(),AliAnalysisManager::kInputContainer);
141
142   // ----- output data -----
143   
144   //slot 0 : default output tree (by default handled by AliAnalysisTaskSE)
145   AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("ctree0", TTree::Class(),AliAnalysisManager::kOutputContainer,"output.root");
146
147   //now comes user's output objects :
148   
149   // output TH1I for event counting
150   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,"output.root");
151   // output Correction Framework Container (for acceptance & efficiency calculations)
152   AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,"output.root");
153   AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("corr0", THnSparseD::Class(),AliAnalysisManager::kOutputContainer,"output.root");
154
155   cinput0->SetData(analysisChain);
156
157   mgr->AddTask(task);
158   mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
159   mgr->ConnectOutput(task,0,coutput0);
160   mgr->ConnectOutput(task,1,coutput1);
161   mgr->ConnectOutput(task,2,coutput2);
162   mgr->ConnectOutput(task,3,coutput3);
163  
164   printf("READY TO RUN\n");
165   //RUN !!!
166   if (mgr->InitAnalysis()) {
167     mgr->PrintStatus();
168     mgr->StartAnalysis("local",analysisChain);
169   }
170
171   benchmark.Stop("AliSingleTrackTask");
172   benchmark.Show("AliSingleTrackTask");
173
174   return kTRUE ;
175 }
176
177 void Load() {
178
179   //load the required aliroot libraries
180   gSystem->Load("libANALYSIS") ;
181   gSystem->Load("libANALYSISalice") ;
182   gSystem->Load("libCORRFW.so") ;
183
184   //compile online the task class
185   gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include -I$ROOTSYS/include");
186   gROOT->LoadMacro("./AliCFTaskForUnfolding.cxx+");
187 }