]>
Commit | Line | Data |
---|---|---|
a5d92cb4 | 1 | //-----------------------------------------------------// |
2 | // // | |
3 | // Base Macro to evaluate pt-eta unfolding // | |
4 | // for charm mesons // | |
5 | // // | |
6 | // Usage: // | |
7 | // // | |
8 | // 1) By default need output.root file from AliCF... // | |
9 | // task // | |
10 | // 2) ./output.root - used to evaluate the eff // | |
11 | // 3) To correct, the data sample should be spitted in // | |
12 | // 2, in this case in GetMeasuredSpectrum() connect // | |
13 | // the second container ./output_data.root // | |
14 | // 4) Set the number of iterations for Baysian unf // | |
15 | // 5) Set min. Chi2 // | |
16 | // 6) Select if unfold at Acceptance level or // | |
17 | // after PPR cuts // | |
18 | // 7) Gaussian error propagation assumed, may be // | |
19 | // errors overestimated - to be studied // | |
20 | // // | |
21 | //-----------------------------------------------------// | |
22 | // // | |
23 | // a.grelli@uu.nl- Utrecht for R.Vernaut example // | |
24 | // // | |
25 | //-----------------------------------------------------// | |
26 | // // | |
27 | // NOTE: by default the macro is computing the // | |
28 | // unfolding after acceptance step. To evaluate // | |
29 | // after PPR switch on the option in the CF task! // | |
30 | // and change in this macro the selection steps // | |
31 | // // | |
32 | //-----------------------------------------------------// | |
33 | ||
34 | #include "TString.h" | |
35 | #include <iostream> | |
36 | #include <fstream> | |
37 | ||
38 | void UnfoldingHF() { | |
39 | ||
40 | // not jet impemented the "systematics" mode! | |
41 | ||
42 | //TString smode(analysis_mode); | |
43 | //TString hist(method); | |
44 | ||
45 | Load(); | |
46 | ||
47 | printf("==================================================================\n"); | |
48 | printf("=========== RUNNING D MESONS UNFOLDING ==========\n"); | |
49 | printf("==================================================================\n"); | |
50 | ||
51 | // --------------------- get variables from container ----------------------// | |
52 | ||
53 | AliCFDataGrid *measured = (AliCFDataGrid*) GetMeasuredSpectrum(); | |
54 | AliCFDataGrid *reconstructed = (AliCFDataGrid*) GetReconstructedSpectrum(); | |
55 | AliCFDataGrid *generated = (AliCFDataGrid*) GetGeneratedSpectrum(); | |
56 | AliCFEffGrid *efficiency = (AliCFEffGrid*) GetEfficiency(); | |
57 | ||
58 | // the response matrix | |
59 | ||
60 | THnSparse *response = (THnSparse*) GetResponseMatrix(); | |
61 | ||
62 | // -------------------------------------------------------------------------// | |
63 | // // | |
64 | // The HF correction framework has 12 dimensions. This is bad for // | |
65 | // pt-eta unfolding so reduce the dimesions from 12 to 2. // | |
66 | // // | |
67 | //--------------------------------------------------------------------------// | |
68 | ||
69 | ||
70 | THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ; | |
71 | THnSparse* data = CreateGuessed(((AliCFGridSparse*)measured->GetData())->GetGrid()); | |
72 | THnSparse* Reco = CreateGuessed(((AliCFGridSparse*)reconstructed->GetData())->GetGrid()) ; | |
73 | ||
74 | ||
75 | // best guess, traslate the 12 dimentions container over 2 (eta and pt) | |
76 | ||
77 | Int_t dim[2]; | |
78 | ||
79 | dim[0] = 0; | |
80 | dim[1] = 1; | |
81 | ||
82 | THnSparse* BestGuess = guessed ->Projection(2,dim); | |
83 | THnSparse* DataSample = data ->Projection(2,dim); | |
84 | THnSparse* Reconstructed = Reco ->Projection(2,dim); | |
85 | THnSparse* SimpleEff = Reco ->Projection(2,dim); | |
86 | ||
87 | SimpleEff->Divide(Reconstructed,BestGuess,1.,1.); | |
88 | ||
89 | ||
90 | // here I do the unfolding and I set the number of iterations and the chi2 | |
91 | ||
92 | AliCFUnfolding unfolding("unfolding","",2,response,SimpleEff,DataSample,BestGuess); | |
93 | unfolding.SetMaxNumberOfIterations(100); | |
94 | unfolding.SetMaxChi2PerDOF(0.00000000005); | |
95 | unfolding.Unfold(); | |
96 | ||
97 | THnSparse* result = unfolding.GetUnfolded(); | |
98 | ||
99 | ||
100 | TH2D* h_gen = generated ->Project(1,0); | |
101 | TH2D* h_meas = DataSample ->Projection(1,0); | |
102 | TH2D* h_reco = Reconstructed ->Projection(1,0); | |
103 | TH2D* h_guessed = BestGuess ->Projection(1,0); | |
104 | TH2D* h_eff = SimpleEff ->Projection(1,0); | |
105 | TH2D* h_unf = result ->Projection(1,0); | |
106 | ||
107 | // Apply simple efficiency correction | |
108 | ||
109 | ||
110 | TH2D* simpleC = (TH2D*)h_meas->Clone(); | |
111 | simpleC->Divide(h_eff); | |
112 | ||
113 | TCanvas * canvas3 = new TCanvas("Unfolded efficiency map","canvas 3",1000,700); | |
114 | ||
115 | canvas3->cd(); | |
116 | ||
117 | TH2D* ratio = (TH2D*)h_unf->Clone(); | |
118 | ratio->SetTitle("corrected using unfolding"); | |
119 | ratio->Divide(h_unf,h_guessed,1,1); | |
120 | ||
121 | ratio->Draw("cont4z"); | |
122 | ||
123 | TCanvas * canvas4 = new TCanvas("Simple Efficiency map","",1000,700); | |
124 | ||
125 | canvas4->cd(); | |
126 | h_meas->SetTitle("simple correction"); | |
127 | h_meas->Divide(h_eff); | |
128 | h_meas->Divide(h_guessed); | |
129 | h_meas->Draw("cont4z"); | |
130 | ||
131 | TCanvas * distribution2 = new TCanvas("dist2","",1000,700); | |
132 | ||
133 | distribution2->cd(); | |
134 | ||
135 | TH1D* gen2 = generated->Project(1); // generated pt | |
136 | gen2->SetLineColor(1); | |
137 | gen2->SetTitle("D0 #eta distribution"); | |
138 | gen2->GetXaxis()->SetTitle("eta"); | |
139 | gen2->SetMarkerStyle(21); | |
140 | gen2->SetMarkerSize(1.); | |
141 | gen2->SetMarkerColor(1); | |
142 | gen2->Draw(); | |
143 | ||
144 | TH1D* meas2 = measured->Project(1); // generated pt | |
145 | meas2->SetTitle("measurements"); | |
146 | meas2->SetLineColor(2); | |
147 | meas2->SetMarkerStyle(22); | |
148 | meas2->SetMarkerSize(1.); | |
149 | meas2->SetMarkerColor(2); | |
150 | meas2->DrawCopy("same"); | |
151 | ||
152 | TH1D* unf2 = result->Projection(1); // generated pt | |
153 | unf2->SetTitle("measurements"); | |
154 | unf2->SetLineColor(3); | |
155 | unf2->SetMarkerStyle(23); | |
156 | unf2->SetMarkerSize(1.); | |
157 | unf2->SetMarkerColor(3); | |
158 | unf2->DrawCopy("SAME"); | |
159 | ||
160 | corr2 = new TH1D(); | |
161 | corr2->Sumw2(); | |
162 | corr2 = simpleC->ProjectionY(); | |
163 | corr2->SetTitle("measurements"); | |
164 | corr2->SetLineColor(4); | |
165 | corr2->SetMarkerStyle(24); | |
166 | corr2->SetMarkerSize(1.); | |
167 | corr2->SetMarkerColor(4); | |
168 | corr2->Sumw2(); | |
169 | corr2->DrawCopy("SAME"); | |
170 | ||
171 | leg = new TLegend(0.1,0.7,0.48,0.9); | |
172 | leg->SetHeader("Distributions"); | |
173 | leg->AddEntry(gen2,"Generated PYTHIA","l"); | |
174 | leg->AddEntry(meas2,"Reconstructed","l"); | |
175 | leg->AddEntry(corr2,"Corrected for eff","l"); | |
176 | leg->AddEntry(unf2,"Unfolded","l"); | |
177 | ||
178 | leg->Draw(); | |
179 | ||
180 | TCanvas * distribution3 = new TCanvas("dist3","",1000,700); | |
181 | ||
182 | distribution3->cd(); | |
183 | ||
184 | ||
185 | TH1D* gen3 = generated->Project(0); // generated pt | |
186 | Int_t Entries = gen3->GetEntries(); | |
187 | gen3->SetLineColor(1); | |
188 | gen3->SetTitle("D0 pt distribution"); | |
189 | gen3->GetXaxis()->SetTitle("pt (GeV/c)"); | |
190 | gen3->SetMarkerStyle(21); | |
191 | gen3->SetMarkerSize(1.); | |
192 | gen3->SetMarkerColor(1); | |
193 | gen3->Draw(); | |
194 | ||
195 | ||
196 | distribution3->cd(); | |
197 | TH1D* meas3 = measured->Project(0); // generated pt | |
198 | meas3->SetTitle("measurements"); | |
199 | meas3->SetLineColor(2); | |
200 | meas3->SetMarkerStyle(22); | |
201 | meas3->SetMarkerSize(1.); | |
202 | meas3->SetMarkerColor(2); | |
203 | meas3->DrawCopy("SAME"); | |
204 | ||
205 | TH1D* unf3 = result->Projection(0); // | |
206 | unf3->SetTitle("unfolded"); | |
207 | unf3->SetLineColor(3); | |
208 | unf3->SetMarkerStyle(23); | |
209 | unf3->SetMarkerSize(1.); | |
210 | unf3->SetMarkerColor(3); | |
211 | unf3->Sumw2(); | |
212 | unf3->DrawCopy("SAME"); | |
213 | ||
214 | TH1D* corr3 = simpleC->ProjectionX(); | |
215 | corr3->SetTitle("corrected"); | |
216 | corr3->SetLineColor(4); | |
217 | corr3->SetMarkerStyle(24); | |
218 | corr3->SetMarkerSize(1.); | |
219 | corr3->SetMarkerColor(4); | |
220 | corr3->Sumw2(); | |
221 | corr3->DrawCopy("SAME"); | |
222 | ||
223 | leg = new TLegend(0.1,0.7,0.48,0.9); | |
224 | leg -> SetHeader("Distributions"); | |
225 | leg -> AddEntry(gen3,"Generated PYTHIA","l"); | |
226 | leg -> AddEntry(meas3,"Reconstructed","l"); | |
227 | leg -> AddEntry(corr3,"Corrected for eff","l"); | |
228 | leg -> AddEntry(unf3,"Unfolded","l"); | |
229 | ||
230 | leg->Draw(); | |
231 | ||
232 | return; | |
233 | } | |
234 | ||
235 | // ==================================================================== | |
236 | ||
237 | ||
238 | void Load(){ | |
239 | gSystem->Load("libANALYSIS"); | |
240 | gSystem->Load("libCORRFW"); | |
241 | gStyle->SetPalette(1); | |
242 | gStyle->SetOptStat(0); | |
243 | } | |
244 | ||
245 | //___________________________________________________________- | |
246 | ||
247 | TObject* GetMeasuredSpectrum() { | |
248 | TFile * f = TFile::Open("output.root","read"); | |
249 | AliCFContainer* c = (AliCFContainer*)f->Get("container"); | |
250 | AliCFDataGrid* data1 = new AliCFDataGrid("data1","",*c); | |
251 | data1->SetMeasured(3); | |
252 | return data1; | |
253 | } | |
254 | TObject* GetReconstructedSpectrum() { | |
255 | TFile * f = TFile::Open("output.root","read"); | |
256 | AliCFContainer* c = (AliCFContainer*)f->Get("container"); | |
257 | AliCFDataGrid* data2 = new AliCFDataGrid("data2","",*c); | |
258 | data2->SetMeasured(3); | |
259 | ||
260 | return data2; | |
261 | } | |
262 | TObject* GetGeneratedSpectrum() { | |
263 | TFile * f = TFile::Open("output.root","read"); | |
264 | AliCFContainer* c = (AliCFContainer*)f->Get("container"); | |
265 | AliCFDataGrid* data3 = new AliCFDataGrid("data3","",*c); | |
266 | data3->SetMeasured(1); | |
267 | return data3; | |
268 | } | |
269 | TObject* GetEfficiency() { | |
270 | TFile * f = TFile::Open("output.root","read"); | |
271 | AliCFContainer* c = (AliCFContainer*)f->Get("container"); | |
272 | AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c); | |
273 | eff->CalculateEfficiency(3,1); | |
274 | return eff; | |
275 | } | |
276 | THnSparse* GetResponseMatrix() { | |
277 | TFile * f = TFile::Open("output.root","read"); | |
278 | THnSparse* h = (THnSparse*)f->Get("correlation"); | |
279 | return h; | |
280 | } | |
281 | ||
282 | THnSparse* GetResponseMatrixPPR() { | |
283 | TFile * f = TFile::Open("output.root","read"); | |
284 | THnSparse* h = (THnSparse*)f->Get("correlationPPR"); | |
285 | return h; | |
286 | } | |
287 | ||
288 | THnSparse* CreateGuessed(const THnSparse* h) { | |
289 | THnSparse* guessed = (THnSparse*) h->Clone(); | |
290 | return guessed ; | |
291 | } | |
292 | ||
293 |