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