0ab29cfa |
1 | /* $Id$ */ |
2 | |
3 | // |
dca331bb |
4 | // script to correct the multiplicity spectrum + helpers |
0ab29cfa |
5 | // |
6 | |
0b4bfd98 |
7 | void SetTPC() |
0a173978 |
8 | { |
9 | gSystem->Load("libPWG0base"); |
0b4bfd98 |
10 | AliMultiplicityCorrection::SetQualityRegions(kFALSE); |
11 | } |
0a173978 |
12 | |
2440928d |
13 | void draw(const char* fileName = "multiplicity.root", const char* folder = "Multiplicity") |
0b4bfd98 |
14 | { |
2440928d |
15 | loadlibs(); |
0a173978 |
16 | |
0b4bfd98 |
17 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder); |
0a173978 |
18 | |
0b4bfd98 |
19 | TFile::Open(fileName); |
20 | mult->LoadHistograms(); |
0a173978 |
21 | mult->DrawHistograms(); |
447c325d |
22 | |
0b4bfd98 |
23 | return; |
24 | |
447c325d |
25 | TH2* hist = (TH2*) gROOT->FindObject("fCorrelation3_zy"); |
26 | canvas = new TCanvas("c1", "c1", 600, 500); |
27 | hist->SetStats(kFALSE); |
28 | hist->Draw("COLZ"); |
29 | hist->SetTitle(";true multiplicity in |#eta| < 2;measured multiplicity in |#eta| < 2"); |
30 | hist->GetYaxis()->SetTitleOffset(1.1); |
31 | gPad->SetRightMargin(0.15); |
32 | gPad->SetLogz(); |
33 | |
34 | canvas->SaveAs("Plot_Correlation.pdf"); |
0ab29cfa |
35 | } |
36 | |
0f67a57c |
37 | void loadlibs() |
9ca6be09 |
38 | { |
0f67a57c |
39 | gSystem->Load("libTree"); |
40 | gSystem->Load("libVMC"); |
41 | |
42 | gSystem->Load("libSTEERBase"); |
43 | gSystem->Load("libANALYSIS"); |
2440928d |
44 | gSystem->Load("libANALYSISalice"); |
9ca6be09 |
45 | gSystem->Load("libPWG0base"); |
0f67a57c |
46 | } |
9ca6be09 |
47 | |
0a0f4adb |
48 | void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityESD.root", Bool_t chi2 = kTRUE, Int_t histID = 2, Bool_t fullPhaseSpace = kFALSE, Float_t beta = 1e3, Int_t eventType = 0 /* AliMultiplicityCorrection::kTrVtx */) |
0f67a57c |
49 | { |
50 | loadlibs(); |
0a0f4adb |
51 | |
0b4bfd98 |
52 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder); |
9ca6be09 |
53 | |
54 | TFile::Open(fileNameMC); |
0b4bfd98 |
55 | mult->LoadHistograms(); |
9ca6be09 |
56 | |
57 | TFile::Open(fileNameESD); |
cfc19dd5 |
58 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); |
447c325d |
59 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID))); |
60 | //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID)); |
9ca6be09 |
61 | |
cfc19dd5 |
62 | mult->SetMultiplicityESD(histID, hist); |
9ca6be09 |
63 | |
447c325d |
64 | // small hack to get around charge conservation for full phase space ;-) |
65 | if (fullPhaseSpace) |
66 | { |
67 | TH1* corr = mult->GetCorrelation(histID + 4); |
dd701109 |
68 | |
447c325d |
69 | for (Int_t i=2; i<=corr->GetNbinsX(); i+=2) |
70 | for (Int_t j=1; j<=corr->GetNbinsY(); ++j) |
71 | { |
72 | corr->SetBinContent(i, j, corr->GetBinContent(i-1, j)); |
73 | corr->SetBinError(i, j, corr->GetBinError(i-1, j)); |
74 | } |
75 | } |
dd701109 |
76 | |
447c325d |
77 | if (chi2) |
78 | { |
0f67a57c |
79 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, beta); |
80 | //mult->SetCreateBigBin(kFALSE); |
81 | //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE); |
82 | //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE); |
83 | //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e5); |
0b4bfd98 |
84 | //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5); |
44466df2 |
85 | //mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kTRUE, hist2->ProjectionY("mymchist")); |
0a0f4adb |
86 | mult->ApplyMinuitFit(histID, fullPhaseSpace, eventType, kFALSE); //hist2->ProjectionY("mymchist")); |
447c325d |
87 | } |
88 | else |
dd701109 |
89 | { |
0a0f4adb |
90 | mult->ApplyBayesianMethod(histID, fullPhaseSpace, eventType, 0.2, 100); |
dd701109 |
91 | } |
92 | |
dca331bb |
93 | TFile* file = TFile::Open("unfolded.root", "RECREATE"); |
144ff489 |
94 | mult->SaveHistograms(); |
95 | file->Write(); |
96 | file->Close(); |
0a0f4adb |
97 | |
98 | mult->DrawComparison((chi2) ? "MinuitChi2" : "Bayesian", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist")); |
99 | } |
100 | |
101 | void CompareChi2Bayesian(Int_t histID = 2, Bool_t showMC = kFALSE, const char* chi2File = "chi2.root", const char* bayesianFile = "bayesian.root", const char* mcFile = "chi2.root") |
102 | { |
103 | const char* folder = "Multiplicity"; |
104 | |
105 | loadlibs(); |
106 | |
107 | AliMultiplicityCorrection* chi2 = new AliMultiplicityCorrection(folder, folder); |
108 | TFile::Open(chi2File); |
109 | chi2->LoadHistograms(); |
110 | |
111 | AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection(folder, folder); |
112 | TFile::Open(bayesianFile); |
113 | bayesian->LoadHistograms(); |
114 | |
115 | histC = chi2->GetMultiplicityESDCorrected(histID); |
116 | histB = bayesian->GetMultiplicityESDCorrected(histID); |
117 | |
118 | c = new TCanvas("CompareChi2Bayesian", "CompareChi2Bayesian", 800, 600); |
119 | c->SetRightMargin(0.05); |
120 | c->SetTopMargin(0.05); |
121 | c->SetLogy(); |
122 | |
123 | histC->SetTitle(";N;P(N)"); |
124 | histC->SetStats(kFALSE); |
125 | histC->GetXaxis()->SetRangeUser(0, 100); |
126 | |
127 | histC->SetLineColor(1); |
128 | histB->SetLineColor(2); |
129 | |
130 | histC->Draw(); |
131 | histB->Draw("SAME"); |
132 | |
133 | legend = new TLegend(0.2, 0.2, 0.4, 0.4); |
134 | legend->SetFillColor(0); |
135 | |
136 | legend->AddEntry(histC, "Chi2"); |
137 | legend->AddEntry(histB, "Bayesian"); |
138 | |
139 | if (showMC) |
140 | { |
141 | AliMultiplicityCorrection* mc = new AliMultiplicityCorrection(folder, folder); |
142 | TFile::Open(mcFile); |
143 | mc->LoadHistograms(); |
144 | |
145 | histMC = mc->GetMultiplicityVtx(histID)->ProjectionY("mc", 1, mc->GetMultiplicityVtx(histID)->GetNbinsX()); |
146 | histMC->Scale(1.0 / histMC->Integral()); |
147 | |
148 | histMC->Draw("SAME"); |
149 | histMC->SetLineColor(4); |
150 | legend->AddEntry(histMC, "MC"); |
151 | } |
152 | |
153 | legend->Draw(); |
154 | |
155 | c->SaveAs(Form("%s.gif", c->GetName())); |
156 | c->SaveAs(Form("%s.eps", c->GetName())); |
447c325d |
157 | } |
158 | |
159 | void* fit2Step(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE) |
160 | { |
161 | gSystem->Load("libPWG0base"); |
162 | |
163 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
164 | |
165 | TFile::Open(fileNameMC); |
166 | mult->LoadHistograms("Multiplicity"); |
167 | |
168 | TFile::Open(fileNameESD); |
169 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); |
170 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID))); |
171 | //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID)); |
172 | |
173 | mult->SetMultiplicityESD(histID, hist); |
174 | |
175 | // small hack to get around charge conservation for full phase space ;-) |
176 | if (fullPhaseSpace) |
177 | { |
178 | TH1* corr = mult->GetCorrelation(histID + 4); |
179 | |
180 | for (Int_t i=2; i<=corr->GetNbinsX(); i+=2) |
181 | for (Int_t j=1; j<=corr->GetNbinsY(); ++j) |
182 | { |
183 | corr->SetBinContent(i, j, corr->GetBinContent(i-1, j)); |
184 | corr->SetBinError(i, j, corr->GetBinError(i-1, j)); |
185 | } |
186 | } |
187 | |
188 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
189 | mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); |
190 | mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist")); |
191 | |
192 | TH1* result = (TH1*) mult->GetMultiplicityESDCorrected((fullPhaseSpace) ? 4 : histID))->Clone("firstresult"); |
193 | |
194 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 100000); |
195 | mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE, result); |
196 | mult->DrawComparison("MinuitChi2_Step2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist")); |
0a173978 |
197 | |
198 | return mult; |
199 | } |
cfc19dd5 |
200 | |
201 | const char* GetRegName(Int_t type) |
202 | { |
203 | switch (type) |
204 | { |
205 | case AliMultiplicityCorrection::kNone: return "None"; break; |
206 | case AliMultiplicityCorrection::kPol0: return "Pol0"; break; |
207 | case AliMultiplicityCorrection::kPol1: return "Pol1"; break; |
208 | case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break; |
209 | case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break; |
0b4bfd98 |
210 | case AliMultiplicityCorrection::kLog : return "Log"; break; |
dd701109 |
211 | } |
212 | return 0; |
213 | } |
214 | |
215 | const char* GetEventTypeName(Int_t type) |
216 | { |
217 | switch (type) |
218 | { |
219 | case AliMultiplicityCorrection::kTrVtx: return "trigger, vertex"; break; |
220 | case AliMultiplicityCorrection::kMB: return "minimum bias"; break; |
221 | case AliMultiplicityCorrection::kINEL: return "inelastic"; break; |
cfc19dd5 |
222 | } |
223 | return 0; |
224 | } |
225 | |
447c325d |
226 | void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3) |
227 | { |
228 | gSystem->mkdir(targetDir); |
229 | |
230 | gSystem->Load("libPWG0base"); |
231 | |
232 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
233 | TFile::Open(fileNameMC); |
234 | mult->LoadHistograms("Multiplicity"); |
235 | |
236 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
237 | TFile::Open(fileNameESD); |
238 | multESD->LoadHistograms("Multiplicity"); |
239 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
240 | |
241 | Int_t count = 0; // just to order the saved images... |
242 | |
243 | TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE"); |
244 | |
0b4bfd98 |
245 | Int_t colors[3] = {1, 2, 4}; |
0f67a57c |
246 | Int_t markers[20] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3, 4, 5, 6}; |
0b4bfd98 |
247 | |
248 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type) |
249 | //for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) |
447c325d |
250 | { |
251 | TString tmp; |
252 | tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type)); |
253 | |
254 | TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600); |
255 | |
0f67a57c |
256 | for (Int_t i = 1; i <= 5; i++) |
447c325d |
257 | { |
0f67a57c |
258 | Int_t iterArray[5] = {5, 20, 50, 100, -1}; |
0b4bfd98 |
259 | //Int_t iter = i * 40 - 20; |
260 | Int_t iter = iterArray[i-1]; |
261 | |
262 | TGraph* fitResultsMC[3]; |
263 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
264 | { |
265 | fitResultsMC[region] = new TGraph; |
266 | fitResultsMC[region]->SetTitle(Form("%d iter. - reg. %d", iter, region+1)); |
267 | fitResultsMC[region]->GetXaxis()->SetTitle("smoothing parameter #alpha"); |
268 | fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region)); |
269 | fitResultsMC[region]->SetName(Form("%s_MC_%d", tmp.Data(), i * AliMultiplicityCorrection::kQualityRegions + region - 2)); |
270 | fitResultsMC[region]->SetFillColor(0); |
0f67a57c |
271 | //fitResultsMC[region]->SetMarkerStyle(markers[(i-1) * AliMultiplicityCorrection::kQualityRegions + region]); |
272 | fitResultsMC[region]->SetMarkerStyle(markers[(i-1)]); |
0b4bfd98 |
273 | fitResultsMC[region]->SetLineColor(colors[region]); |
274 | } |
275 | |
447c325d |
276 | TGraph* fitResultsRes = new TGraph; |
277 | fitResultsRes->SetTitle(Form("%d iterations", iter)); |
278 | fitResultsRes->SetName(Form("%s_Res_%d", tmp.Data(), i)); |
279 | fitResultsRes->GetXaxis()->SetTitle("smoothing parameter"); |
280 | fitResultsRes->GetYaxis()->SetTitle("P_{2}"); |
281 | |
447c325d |
282 | fitResultsRes->SetFillColor(0); |
447c325d |
283 | fitResultsRes->SetMarkerStyle(19+i); |
0b4bfd98 |
284 | fitResultsRes->SetMarkerColor(1); |
285 | fitResultsRes->SetLineColor(1); |
447c325d |
286 | |
287 | for (Float_t weight = 0.0; weight < 1.01; weight += 0.2) |
288 | { |
289 | Float_t chi2MC = 0; |
290 | Float_t residuals = 0; |
291 | |
0b4bfd98 |
292 | mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter, 0, kFALSE); |
447c325d |
293 | mult->DrawComparison(Form("%s/BayesianIterSmooth_%03d_%d_%d_%f", targetDir, count++, type, iter, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); |
294 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
295 | |
0b4bfd98 |
296 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
297 | fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region)); |
298 | |
447c325d |
299 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); |
300 | } |
301 | |
0b4bfd98 |
302 | graphFile->cd(); |
303 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
304 | fitResultsMC[region]->Write(); |
305 | |
447c325d |
306 | fitResultsRes->Write(); |
307 | } |
308 | } |
0b4bfd98 |
309 | |
310 | graphFile->Close(); |
447c325d |
311 | } |
312 | |
313 | void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE) |
314 | { |
315 | gSystem->Load("libPWG0base"); |
316 | |
317 | TString name; |
318 | TFile* graphFile = 0; |
319 | if (type == -1) |
320 | { |
321 | name = "EvaluateChi2Method"; |
322 | graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir)); |
323 | } |
324 | else |
325 | { |
326 | name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type)); |
327 | graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir)); |
328 | } |
329 | |
330 | TCanvas* canvas = new TCanvas(name, name, 800, 500); |
331 | if (type == -1) |
0b4bfd98 |
332 | { |
447c325d |
333 | canvas->SetLogx(); |
0b4bfd98 |
334 | canvas->SetLogy(); |
335 | } |
336 | canvas->SetTopMargin(0.05); |
337 | canvas->SetGridx(); |
338 | canvas->SetGridy(); |
447c325d |
339 | |
0b4bfd98 |
340 | TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98); |
447c325d |
341 | legend->SetFillColor(0); |
342 | |
343 | Int_t count = 1; |
344 | |
345 | Float_t xMin = 1e20; |
346 | Float_t xMax = 0; |
347 | |
348 | Float_t yMin = 1e20; |
349 | Float_t yMax = 0; |
350 | |
0b4bfd98 |
351 | Float_t yMinRegion[3]; |
352 | for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i) |
353 | yMinRegion[i] = 1e20; |
354 | |
447c325d |
355 | TString xaxis, yaxis; |
356 | |
357 | while (1) |
358 | { |
359 | TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count)); |
360 | TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count)); |
361 | |
0b4bfd98 |
362 | if (!mc) |
447c325d |
363 | break; |
364 | |
365 | xaxis = mc->GetXaxis()->GetTitle(); |
366 | yaxis = mc->GetYaxis()->GetTitle(); |
367 | |
368 | mc->Print(); |
447c325d |
369 | |
0b4bfd98 |
370 | if (res) |
371 | res->Print(); |
447c325d |
372 | |
0b4bfd98 |
373 | xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin()); |
374 | yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin()); |
447c325d |
375 | |
0b4bfd98 |
376 | xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax()); |
377 | yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax()); |
378 | |
379 | if (plotRes && res) |
447c325d |
380 | { |
0b4bfd98 |
381 | xMin = TMath::Min(xMin, res->GetXaxis()->GetXmin()); |
382 | yMin = TMath::Min(yMin, res->GetYaxis()->GetXmin()); |
447c325d |
383 | |
0b4bfd98 |
384 | xMax = TMath::Max(xMax, res->GetXaxis()->GetXmax()); |
385 | yMax = TMath::Max(yMax, res->GetYaxis()->GetXmax()); |
447c325d |
386 | } |
0b4bfd98 |
387 | |
388 | for (Int_t i=0; i<mc->GetN(); ++i) |
389 | yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]); |
390 | |
391 | count++; |
447c325d |
392 | } |
393 | |
0b4bfd98 |
394 | for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i) |
395 | Printf("Minimum for region %d is %f", i, yMinRegion[i]); |
396 | |
447c325d |
397 | if (type >= 0) |
398 | { |
399 | xaxis = "smoothing parameter"; |
400 | } |
401 | else if (type == -1) |
402 | { |
403 | xaxis = "weight parameter"; |
0b4bfd98 |
404 | xMax *= 5; |
447c325d |
405 | } |
0b4bfd98 |
406 | //yaxis = "P_{1} (2 <= t <= 150)"; |
447c325d |
407 | |
408 | printf("%f %f %f %f\n", xMin, xMax, yMin, yMax); |
409 | |
410 | TGraph* dummy = new TGraph; |
411 | dummy->SetPoint(0, xMin, yMin); |
412 | dummy->SetPoint(1, xMax, yMax); |
413 | dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data())); |
414 | |
415 | dummy->SetMarkerColor(0); |
416 | dummy->Draw("AP"); |
0b4bfd98 |
417 | dummy->GetYaxis()->SetMoreLogLabels(1); |
447c325d |
418 | |
419 | count = 1; |
420 | |
421 | while (1) |
422 | { |
423 | TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count)); |
424 | TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count)); |
425 | |
0b4bfd98 |
426 | //printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res); |
447c325d |
427 | |
0b4bfd98 |
428 | if (!mc) |
447c325d |
429 | break; |
430 | |
431 | printf("Loaded %d sucessful.\n", count); |
432 | |
44466df2 |
433 | if (type == -1) |
447c325d |
434 | { |
44466df2 |
435 | legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3)); |
447c325d |
436 | } |
44466df2 |
437 | else |
447c325d |
438 | legend->AddEntry(mc); |
439 | |
440 | mc->Draw("SAME PC"); |
441 | |
0b4bfd98 |
442 | if (plotRes && res) |
447c325d |
443 | { |
444 | legend->AddEntry(res); |
445 | res->Draw("SAME PC"); |
446 | } |
447 | |
448 | count++; |
449 | } |
450 | |
451 | legend->Draw(); |
452 | |
453 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
454 | canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName())); |
455 | } |
456 | |
0f67a57c |
457 | void EvaluateDrawResultRegions(const char* targetDir, Int_t type = 0) |
458 | { |
459 | gSystem->Load("libPWG0base"); |
460 | |
461 | TString name; |
462 | TFile* graphFile = 0; |
463 | if (type == -1) |
464 | { |
465 | name = "EvaluateChi2Method"; |
466 | graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir)); |
467 | } |
468 | else |
469 | { |
470 | name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type)); |
471 | graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir)); |
472 | } |
473 | |
474 | TCanvas* canvas = new TCanvas(name, name, 800, 1200); |
475 | //canvas->Divide(1, AliMultiplicityCorrection::kQualityRegions, 0, 0); |
476 | canvas->Range(0, 0, 1, 1); |
477 | |
478 | TPad* pad[3]; |
479 | pad[0] = new TPad(Form("%s_pad1", name.Data()), "", 0.02, 0.05, 0.98, 0.35); |
480 | pad[1] = new TPad(Form("%s_pad2", name.Data()), "", 0.02, 0.35, 0.98, 0.65); |
481 | pad[2] = new TPad(Form("%s_pad3", name.Data()), "", 0.02, 0.65, 0.98, 0.95); |
482 | |
483 | Float_t yMinRegion[3]; |
484 | for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i) |
485 | yMinRegion[i] = 1e20; |
486 | |
487 | for (Int_t region = 1; region <= AliMultiplicityCorrection::kQualityRegions; region++) |
488 | { |
489 | canvas->cd(); |
490 | pad[region-1]->Draw(); |
491 | pad[region-1]->cd(); |
492 | pad[region-1]->SetRightMargin(0.05); |
493 | |
494 | if (region != 1) |
495 | pad[region-1]->SetBottomMargin(0); |
496 | if (region != AliMultiplicityCorrection::kQualityRegions) |
497 | pad[region-1]->SetTopMargin(0); |
498 | |
499 | if (type == -1) |
500 | { |
501 | pad[region-1]->SetLogx(); |
502 | pad[region-1]->SetLogy(); |
503 | } |
504 | pad[region-1]->SetTopMargin(0.05); |
505 | pad[region-1]->SetGridx(); |
506 | pad[region-1]->SetGridy(); |
507 | |
508 | TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98); |
509 | legend->SetFillColor(0); |
510 | |
511 | Int_t count = 0; |
512 | |
513 | Float_t xMin = 1e20; |
514 | Float_t xMax = 0; |
515 | |
516 | Float_t yMin = 1e20; |
517 | Float_t yMax = 0; |
518 | |
519 | TString xaxis, yaxis; |
520 | |
521 | while (1) |
522 | { |
523 | count++; |
524 | |
525 | TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count)); |
526 | if (!mc) |
527 | break; |
528 | |
529 | if (TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE) |
530 | continue; |
531 | |
532 | xaxis = mc->GetXaxis()->GetTitle(); |
533 | yaxis = mc->GetYaxis()->GetTitle(); |
534 | |
535 | mc->Print(); |
536 | |
537 | xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin()); |
538 | yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin()); |
539 | |
540 | xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax()); |
541 | yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax()); |
542 | |
543 | for (Int_t i=0; i<mc->GetN(); ++i) |
544 | yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]); |
545 | } |
546 | |
547 | if (type >= 0) |
548 | { |
549 | xaxis = "smoothing parameter"; |
550 | } |
551 | else if (type == -1) |
552 | { |
553 | xaxis = "weight parameter"; |
554 | xMax *= 5; |
555 | } |
556 | yaxis = "P_{1}"; // (2 <= t <= 150)"; |
557 | |
558 | printf("%f %f %f %f\n", xMin, xMax, yMin, yMax); |
559 | |
560 | TGraph* dummy = new TGraph; |
561 | dummy->SetPoint(0, xMin, yMin); |
562 | dummy->SetPoint(1, xMax, yMax); |
563 | dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data())); |
564 | |
565 | dummy->SetMarkerColor(0); |
566 | dummy->Draw("AP"); |
567 | dummy->GetYaxis()->SetMoreLogLabels(1); |
568 | |
569 | count = 0; |
570 | |
571 | while (1) |
572 | { |
573 | count++; |
574 | |
575 | TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count)); |
576 | if (!mc) |
577 | break; |
578 | if (TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE) |
579 | continue; |
580 | |
581 | printf("Loaded %d sucessful.\n", count); |
582 | |
583 | if (type == -1) |
584 | { |
585 | legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3)); |
586 | } |
587 | else |
588 | legend->AddEntry(mc); |
589 | |
590 | mc->Draw("SAME PC"); |
591 | |
592 | } |
593 | |
594 | legend->Draw(); |
595 | } |
596 | |
597 | for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i) |
598 | Printf("Minimum for region %d is %f", i, yMinRegion[i]); |
599 | |
600 | canvas->Modified(); |
601 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
602 | canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName())); |
603 | } |
604 | |
0b4bfd98 |
605 | void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", const char* targetDir, Int_t histID = 3) |
dd701109 |
606 | { |
607 | gSystem->mkdir(targetDir); |
608 | |
cfc19dd5 |
609 | gSystem->Load("libPWG0base"); |
610 | |
611 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
612 | |
613 | TFile::Open(fileNameMC); |
614 | mult->LoadHistograms("Multiplicity"); |
615 | |
616 | TFile::Open(fileNameESD); |
617 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); |
618 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID)); |
619 | |
620 | mult->SetMultiplicityESD(histID, hist); |
621 | |
dd701109 |
622 | Int_t count = 0; // just to order the saved images... |
0b4bfd98 |
623 | Int_t colors[3] = {1, 2, 4}; |
624 | Int_t markers[12] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3}; |
dd701109 |
625 | |
447c325d |
626 | TGraph* fitResultsRes = 0; |
627 | |
628 | TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE"); |
cfc19dd5 |
629 | |
0f67a57c |
630 | for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type) |
0b4bfd98 |
631 | // for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kEntropy; type <= AliMultiplicityCorrection::kEntropy; ++type) |
0f67a57c |
632 | //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type) |
cfc19dd5 |
633 | { |
0b4bfd98 |
634 | TGraph* fitResultsMC[3]; |
635 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
636 | { |
637 | fitResultsMC[region] = new TGraph; |
44466df2 |
638 | fitResultsMC[region]->SetTitle(Form("Eq. (%d) - reg. %d", type+9, region+1)); |
0b4bfd98 |
639 | fitResultsMC[region]->GetXaxis()->SetTitle("weight parameter #alpha"); |
640 | fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region)); |
641 | fitResultsMC[region]->SetName(Form("EvaluateChi2Method_MC_%d", type * AliMultiplicityCorrection::kQualityRegions + region - 2)); |
642 | fitResultsMC[region]->SetFillColor(0); |
643 | fitResultsMC[region]->SetMarkerStyle(markers[(type-1) * AliMultiplicityCorrection::kQualityRegions + region]); |
644 | fitResultsMC[region]->SetLineColor(colors[region]); |
645 | } |
646 | |
447c325d |
647 | fitResultsRes = new TGraph; |
648 | fitResultsRes->SetTitle(Form("%s residual chi2", GetRegName(type))); |
649 | fitResultsRes->SetName(Form("EvaluateChi2Method_Res_%d", type)); |
650 | fitResultsRes->GetXaxis()->SetTitle("Weight Parameter"); |
cfc19dd5 |
651 | |
cfc19dd5 |
652 | fitResultsRes->SetFillColor(0); |
dd701109 |
653 | fitResultsRes->SetMarkerStyle(23+type); |
cfc19dd5 |
654 | fitResultsRes->SetMarkerColor(kRed); |
655 | fitResultsRes->SetLineColor(kRed); |
656 | |
0b4bfd98 |
657 | for (Int_t i=0; i<7; ++i) |
447c325d |
658 | { |
659 | Float_t weight = TMath::Power(TMath::Sqrt(10), i+6); |
660 | //Float_t weight = TMath::Power(10, i+2); |
cfc19dd5 |
661 | |
447c325d |
662 | //if (type == AliMultiplicityCorrection::kEntropy) weight = 1e4 * (i+1) * 1.5; |
cfc19dd5 |
663 | |
cfc19dd5 |
664 | Float_t chi2MC = 0; |
665 | Float_t residuals = 0; |
447c325d |
666 | Float_t chi2Limit = 0; |
667 | |
668 | TString runName; |
669 | runName.Form("MinuitChi2_%02d_%d_%f", count++, type, weight); |
cfc19dd5 |
670 | |
671 | mult->SetRegularizationParameters(type, weight); |
447c325d |
672 | Int_t status = mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
673 | mult->DrawComparison(Form("%s/%s", targetDir, runName.Data()), histID, kFALSE, kTRUE, hist2->ProjectionY()); |
674 | if (status != 0) |
675 | { |
676 | printf("MINUIT did not succeed. Skipping...\n"); |
677 | continue; |
678 | } |
679 | |
dd701109 |
680 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
447c325d |
681 | TH1* result = mult->GetMultiplicityESDCorrected(histID); |
682 | result->SetName(runName); |
683 | result->Write(); |
cfc19dd5 |
684 | |
0b4bfd98 |
685 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
686 | fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region)); |
687 | |
cfc19dd5 |
688 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); |
cfc19dd5 |
689 | } |
690 | |
447c325d |
691 | graphFile->cd(); |
0b4bfd98 |
692 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
693 | fitResultsMC[region]->Write(); |
447c325d |
694 | fitResultsRes->Write(); |
cfc19dd5 |
695 | } |
696 | |
447c325d |
697 | graphFile->Close(); |
dd701109 |
698 | } |
699 | |
700 | void EvaluateChi2MethodAll() |
701 | { |
702 | EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); |
703 | EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); |
704 | EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); |
705 | EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); |
706 | EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); |
707 | } |
708 | |
709 | void EvaluateBayesianMethodAll() |
710 | { |
711 | EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); |
712 | EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); |
713 | EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); |
714 | EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); |
715 | EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); |
716 | } |
717 | |
447c325d |
718 | void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3) |
dd701109 |
719 | { |
720 | gSystem->mkdir(targetDir); |
721 | |
722 | gSystem->Load("libPWG0base"); |
723 | |
724 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
725 | |
726 | TFile::Open(fileNameMC); |
727 | mult->LoadHistograms("Multiplicity"); |
728 | |
729 | TFile::Open(fileNameESD); |
730 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
731 | multESD->LoadHistograms("Multiplicity"); |
732 | |
733 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
734 | |
447c325d |
735 | TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 1200); |
736 | canvas->Divide(3, 3); |
dd701109 |
737 | |
738 | Int_t count = 0; |
739 | |
447c325d |
740 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type) |
dd701109 |
741 | { |
447c325d |
742 | TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY("mymc"); |
743 | mc->Sumw2(); |
dd701109 |
744 | |
447c325d |
745 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
dd701109 |
746 | mult->ApplyMinuitFit(histID, kFALSE, type); |
747 | mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc); |
748 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result"); |
749 | |
750 | mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1); |
751 | mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc); |
752 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult"); |
753 | |
754 | mc->GetXaxis()->SetRangeUser(0, 150); |
755 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
756 | |
447c325d |
757 | /* // skip errors for now |
dd701109 |
758 | for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i) |
759 | { |
760 | chi2Result->SetBinError(i, 0); |
761 | bayesResult->SetBinError(i, 0); |
447c325d |
762 | }*/ |
dd701109 |
763 | |
764 | canvas->cd(++count); |
765 | mc->SetFillColor(kYellow); |
766 | mc->DrawCopy(); |
767 | chi2Result->SetLineColor(kRed); |
768 | chi2Result->DrawCopy("SAME"); |
769 | bayesResult->SetLineColor(kBlue); |
770 | bayesResult->DrawCopy("SAME"); |
771 | gPad->SetLogy(); |
772 | |
773 | canvas->cd(count + 3); |
447c325d |
774 | chi2ResultRatio = (TH1*) chi2Result->Clone("chi2ResultRatio"); |
775 | bayesResultRatio = (TH1*) bayesResult->Clone("bayesResultRatio"); |
776 | chi2ResultRatio->Divide(chi2Result, mc, 1, 1, ""); |
777 | bayesResultRatio->Divide(bayesResult, mc, 1, 1, ""); |
dd701109 |
778 | |
447c325d |
779 | chi2ResultRatio->GetYaxis()->SetRangeUser(0.5, 1.5); |
dd701109 |
780 | |
447c325d |
781 | chi2ResultRatio->DrawCopy("HIST"); |
782 | bayesResultRatio->DrawCopy("SAME HIST"); |
dd701109 |
783 | |
447c325d |
784 | canvas->cd(count + 6); |
785 | chi2Result->Divide(chi2Result, bayesResult, 1, 1, ""); |
786 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
787 | chi2Result->DrawCopy("HIST"); |
dd701109 |
788 | } |
789 | |
790 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
791 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); |
792 | } |
793 | |
447c325d |
794 | void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3) |
dd701109 |
795 | { |
796 | gSystem->Load("libPWG0base"); |
797 | |
798 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
799 | |
800 | TFile::Open(fileNameMC); |
801 | mult->LoadHistograms("Multiplicity"); |
802 | |
447c325d |
803 | const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" }; |
dd701109 |
804 | |
0b4bfd98 |
805 | TGraph* fitResultsChi2[3]; |
806 | TGraph* fitResultsBayes[3]; |
807 | |
808 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
809 | { |
810 | fitResultsChi2[region] = new TGraph; |
811 | fitResultsChi2[region]->SetTitle(";Nevents;Chi2"); |
812 | fitResultsChi2[region]->SetName(Form("fitResultsChi2_%d", region)); |
813 | fitResultsChi2[region]->SetMarkerStyle(20+region); |
814 | |
815 | fitResultsBayes[region] = new TGraph; |
816 | fitResultsBayes[region]->SetTitle(";Nevents;Chi2"); |
817 | fitResultsBayes[region]->SetName(Form("fitResultsBayes_%d", region)); |
818 | fitResultsBayes[region]->SetMarkerStyle(20+region); |
819 | fitResultsBayes[region]->SetMarkerColor(2); |
820 | } |
821 | |
447c325d |
822 | TGraph* fitResultsChi2Limit = new TGraph; fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach"); |
823 | TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach"); |
0b4bfd98 |
824 | TGraph* fitResultsChi2Res = new TGraph; fitResultsChi2Res->SetTitle(";Nevents;Chi2"); |
825 | TGraph* fitResultsBayesRes = new TGraph; fitResultsBayesRes->SetTitle(";Nevents;Chi2"); |
447c325d |
826 | |
447c325d |
827 | fitResultsChi2Limit->SetName("fitResultsChi2Limit"); |
828 | fitResultsBayesLimit->SetName("fitResultsBayesLimit"); |
0b4bfd98 |
829 | fitResultsChi2Res->SetName("fitResultsChi2Res"); |
830 | fitResultsBayesRes->SetName("fitResultsBayesRes"); |
dd701109 |
831 | |
832 | TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600); |
833 | canvas->Divide(5, 2); |
834 | |
835 | Float_t min = 1e10; |
836 | Float_t max = 0; |
837 | |
447c325d |
838 | TFile* file = TFile::Open("StatisticsPlot.root", "RECREATE"); |
839 | file->Close(); |
840 | |
dd701109 |
841 | for (Int_t i=0; i<5; ++i) |
842 | { |
843 | TFile::Open(files[i]); |
844 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
845 | multESD->LoadHistograms("Multiplicity"); |
846 | |
847 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
848 | Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries(); |
447c325d |
849 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i)); |
dd701109 |
850 | |
447c325d |
851 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
dd701109 |
852 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
853 | mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); |
854 | |
dd701109 |
855 | Int_t chi2MCLimit = 0; |
0b4bfd98 |
856 | Float_t chi2Residuals = 0; |
857 | mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals); |
858 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
859 | { |
860 | fitResultsChi2[region]->SetPoint(fitResultsChi2[region]->GetN(), nEntries, mult->GetQuality(region)); |
861 | min = TMath::Min(min, mult->GetQuality(region)); |
862 | max = TMath::Max(max, mult->GetQuality(region)); |
863 | } |
dd701109 |
864 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit); |
0b4bfd98 |
865 | fitResultsChi2Res->SetPoint(fitResultsChi2Res->GetN(), nEntries, chi2Residuals); |
dd701109 |
866 | |
447c325d |
867 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i)); |
dd701109 |
868 | |
0b4bfd98 |
869 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0, kFALSE); |
dd701109 |
870 | mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
447c325d |
871 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i)); |
0b4bfd98 |
872 | mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals); |
873 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
874 | { |
875 | fitResultsBayes[region]->SetPoint(fitResultsBayes[region]->GetN(), nEntries, mult->GetQuality(region)); |
876 | min = TMath::Min(min, mult->GetQuality(region)); |
877 | max = TMath::Max(max, mult->GetQuality(region)); |
878 | } |
dd701109 |
879 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit); |
0b4bfd98 |
880 | fitResultsBayesRes->SetPoint(fitResultsBayesRes->GetN(), nEntries, chi2Residuals); |
dd701109 |
881 | |
dd701109 |
882 | mc->GetXaxis()->SetRangeUser(0, 150); |
883 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
884 | |
885 | // skip errors for now |
886 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
887 | { |
888 | chi2Result->SetBinError(j, 0); |
889 | bayesResult->SetBinError(j, 0); |
890 | } |
891 | |
892 | canvas->cd(i+1); |
893 | mc->SetFillColor(kYellow); |
894 | mc->DrawCopy(); |
895 | chi2Result->SetLineColor(kRed); |
896 | chi2Result->DrawCopy("SAME"); |
897 | bayesResult->SetLineColor(kBlue); |
898 | bayesResult->DrawCopy("SAME"); |
899 | gPad->SetLogy(); |
900 | |
901 | canvas->cd(i+6); |
902 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); |
903 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); |
904 | |
905 | // skip errors for now |
906 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
907 | { |
908 | chi2Result->SetBinError(j, 0); |
909 | bayesResult->SetBinError(j, 0); |
910 | } |
911 | |
912 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); |
913 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
914 | |
915 | chi2Result->DrawCopy(""); |
916 | bayesResult->DrawCopy("SAME"); |
447c325d |
917 | |
918 | TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE"); |
919 | mc->Write(); |
920 | chi2Result->Write(); |
921 | bayesResult->Write(); |
922 | file->Close(); |
dd701109 |
923 | } |
924 | |
925 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); |
926 | canvas->SaveAs(Form("%s.C", canvas->GetName())); |
927 | |
928 | TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400); |
929 | canvas2->Divide(2, 1); |
930 | |
931 | canvas2->cd(1); |
dd701109 |
932 | |
0b4bfd98 |
933 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
934 | { |
935 | fitResultsChi2[region]->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); |
936 | fitResultsChi2[region]->Draw(((region == 0) ? "AP" : "P SAME")); |
937 | |
938 | fitResultsBayes[region]->Draw("P SAME"); |
939 | } |
dd701109 |
940 | |
941 | gPad->SetLogy(); |
942 | |
943 | canvas2->cd(2); |
944 | fitResultsChi2Limit->SetMarkerStyle(20); |
945 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); |
946 | fitResultsChi2Limit->Draw("AP"); |
947 | |
948 | fitResultsBayesLimit->SetMarkerStyle(3); |
949 | fitResultsBayesLimit->SetMarkerColor(2); |
950 | fitResultsBayesLimit->Draw("P SAME"); |
951 | |
952 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); |
953 | canvas2->SaveAs(Form("%s.C", canvas2->GetName())); |
447c325d |
954 | |
955 | TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE"); |
0b4bfd98 |
956 | |
957 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
958 | { |
959 | fitResultsChi2[region]->Write(); |
960 | fitResultsBayes[region]->Write(); |
961 | } |
447c325d |
962 | fitResultsChi2Limit->Write(); |
963 | fitResultsBayesLimit->Write(); |
0b4bfd98 |
964 | fitResultsChi2Res->Write(); |
965 | fitResultsBayesRes->Write(); |
447c325d |
966 | file->Close(); |
dd701109 |
967 | } |
968 | |
447c325d |
969 | void StartingConditions(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3) |
dd701109 |
970 | { |
971 | gSystem->Load("libPWG0base"); |
972 | |
973 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
974 | |
975 | TFile::Open(fileNameMC); |
976 | mult->LoadHistograms("Multiplicity"); |
977 | |
447c325d |
978 | const char* files[] = { "multiplicityMC_1M_3.root", "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root" } |
dd701109 |
979 | |
980 | // this one we try to unfold |
981 | TFile::Open(files[0]); |
982 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
983 | multESD->LoadHistograms("Multiplicity"); |
984 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
447c325d |
985 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY("mc"); |
dd701109 |
986 | |
987 | TGraph* fitResultsChi2 = new TGraph; |
988 | fitResultsChi2->SetTitle(";Input Dist ID;Chi2"); |
989 | TGraph* fitResultsBayes = new TGraph; |
990 | fitResultsBayes->SetTitle(";Input Dist ID;Chi2"); |
991 | TGraph* fitResultsChi2Limit = new TGraph; |
992 | fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach"); |
993 | TGraph* fitResultsBayesLimit = new TGraph; |
994 | fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach"); |
995 | |
996 | TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600); |
997 | canvas->Divide(8, 2); |
998 | |
999 | TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400); |
1000 | canvas3->Divide(2, 1); |
1001 | |
1002 | Float_t min = 1e10; |
1003 | Float_t max = 0; |
1004 | |
1005 | TH1* firstChi = 0; |
1006 | TH1* firstBayesian = 0; |
1007 | TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond"); |
1008 | |
1009 | TLegend* legend = new TLegend(0.7, 0.7, 1, 1); |
1010 | |
447c325d |
1011 | TFile* file = TFile::Open("StartingConditions.root", "RECREATE"); |
1012 | mc->Write(); |
1013 | file->Close(); |
1014 | |
dd701109 |
1015 | for (Int_t i=0; i<8; ++i) |
1016 | { |
1017 | if (i == 0) |
1018 | { |
1019 | startCond = (TH1*) mc->Clone("startCond2"); |
1020 | } |
1021 | else if (i < 6) |
1022 | { |
1023 | TFile::Open(files[i-1]); |
1024 | AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2"); |
1025 | multESD2->LoadHistograms("Multiplicity"); |
1026 | startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond"); |
1027 | } |
1028 | else if (i == 6) |
1029 | { |
1030 | func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50); |
1031 | func->SetParNames("scaling", "averagen", "k"); |
1032 | func->SetParLimits(0, 1e-3, 1e10); |
1033 | func->SetParLimits(1, 0.001, 1000); |
1034 | func->SetParLimits(2, 0.001, 1000); |
1035 | |
1036 | func->SetParameters(1, 10, 2); |
1037 | for (Int_t j=2; j<=startCond->GetNbinsX(); j++) |
1038 | startCond->SetBinContent(j, func->Eval(j-1)); |
1039 | } |
1040 | else if (i == 7) |
1041 | { |
1042 | for (Int_t j=1; j<=startCond->GetNbinsX(); j++) |
1043 | startCond->SetBinContent(j, 1); |
1044 | } |
1045 | |
447c325d |
1046 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
dd701109 |
1047 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond); |
1048 | mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); |
1049 | |
1050 | Float_t chi2MC = 0; |
1051 | Int_t chi2MCLimit = 0; |
1052 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
1053 | fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC); |
1054 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit); |
1055 | min = TMath::Min(min, chi2MC); |
1056 | max = TMath::Max(max, chi2MC); |
1057 | |
447c325d |
1058 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i)); |
dd701109 |
1059 | if (!firstChi) |
1060 | firstChi = (TH1*) chi2Result->Clone("firstChi"); |
1061 | |
447c325d |
1062 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, startCond); |
dd701109 |
1063 | mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
447c325d |
1064 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i)); |
dd701109 |
1065 | if (!firstBayesian) |
1066 | firstBayesian = (TH1*) bayesResult->Clone("firstBayesian"); |
1067 | |
1068 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
1069 | fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC); |
1070 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit); |
1071 | |
447c325d |
1072 | TFile* file = TFile::Open("StartingConditions.root", "UPDATE"); |
1073 | chi2Result->Write(); |
1074 | bayesResult->Write(); |
1075 | file->Close(); |
1076 | |
dd701109 |
1077 | min = TMath::Min(min, chi2MC); |
1078 | max = TMath::Max(max, chi2MC); |
1079 | mc->GetXaxis()->SetRangeUser(0, 150); |
1080 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
1081 | |
1082 | // skip errors for now |
1083 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
1084 | { |
1085 | chi2Result->SetBinError(j, 0); |
1086 | bayesResult->SetBinError(j, 0); |
1087 | } |
1088 | |
1089 | canvas3->cd(1); |
1090 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); |
1091 | tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); |
1092 | tmp->Divide(firstChi); |
1093 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1094 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1095 | tmp->SetLineColor(i+1); |
1096 | legend->AddEntry(tmp, Form("%d", i)); |
1097 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1098 | |
1099 | canvas3->cd(2); |
1100 | tmp = (TH1*) bayesResult->Clone("tmp"); |
1101 | tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); |
1102 | tmp->Divide(firstBayesian); |
1103 | tmp->SetLineColor(i+1); |
1104 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1105 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1106 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1107 | |
1108 | canvas->cd(i+1); |
1109 | mc->SetFillColor(kYellow); |
1110 | mc->DrawCopy(); |
1111 | chi2Result->SetLineColor(kRed); |
1112 | chi2Result->DrawCopy("SAME"); |
1113 | bayesResult->SetLineColor(kBlue); |
1114 | bayesResult->DrawCopy("SAME"); |
1115 | gPad->SetLogy(); |
1116 | |
1117 | canvas->cd(i+9); |
1118 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); |
1119 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); |
1120 | |
1121 | // skip errors for now |
1122 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
1123 | { |
1124 | chi2Result->SetBinError(j, 0); |
1125 | bayesResult->SetBinError(j, 0); |
1126 | } |
1127 | |
1128 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); |
1129 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
1130 | |
1131 | chi2Result->DrawCopy(""); |
1132 | bayesResult->DrawCopy("SAME"); |
1133 | } |
1134 | |
1135 | canvas3->cd(1); |
1136 | legend->Draw(); |
1137 | |
cfc19dd5 |
1138 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); |
dd701109 |
1139 | |
1140 | TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400); |
1141 | canvas2->Divide(2, 1); |
1142 | |
1143 | canvas2->cd(1); |
1144 | fitResultsChi2->SetMarkerStyle(20); |
1145 | fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); |
1146 | fitResultsChi2->Draw("AP"); |
1147 | |
1148 | fitResultsBayes->SetMarkerStyle(3); |
1149 | fitResultsBayes->SetMarkerColor(2); |
1150 | fitResultsBayes->Draw("P SAME"); |
1151 | |
1152 | gPad->SetLogy(); |
1153 | |
1154 | canvas2->cd(2); |
1155 | fitResultsChi2Limit->SetMarkerStyle(20); |
1156 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); |
1157 | fitResultsChi2Limit->Draw("AP"); |
1158 | |
1159 | fitResultsBayesLimit->SetMarkerStyle(3); |
1160 | fitResultsBayesLimit->SetMarkerColor(2); |
1161 | fitResultsBayesLimit->Draw("P SAME"); |
1162 | |
1163 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); |
1164 | canvas3->SaveAs(Form("%s.gif", canvas3->GetName())); |
cfc19dd5 |
1165 | } |
1166 | |
447c325d |
1167 | void DifferentSamples(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3) |
1168 | { |
1169 | gSystem->Load("libPWG0base"); |
1170 | |
1171 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1172 | |
1173 | TFile::Open(fileNameMC); |
1174 | mult->LoadHistograms("Multiplicity"); |
1175 | |
1176 | const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" }; |
1177 | |
1178 | TGraph* fitResultsChi2 = new TGraph; |
1179 | fitResultsChi2->SetTitle(";Input Dist ID;Chi2"); |
1180 | TGraph* fitResultsBayes = new TGraph; |
1181 | fitResultsBayes->SetTitle(";Input Dist ID;Chi2"); |
1182 | TGraph* fitResultsChi2Limit = new TGraph; |
1183 | fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach"); |
1184 | TGraph* fitResultsBayesLimit = new TGraph; |
1185 | fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach"); |
1186 | |
1187 | TCanvas* canvasA = new TCanvas("DifferentSamplesA", "DifferentSamplesA", 1200, 600); |
1188 | canvasA->Divide(4, 2); |
1189 | |
1190 | TCanvas* canvasB = new TCanvas("DifferentSamplesB", "DifferentSamplesB", 1200, 600); |
1191 | canvasB->Divide(4, 2); |
1192 | |
1193 | TCanvas* canvas4 = new TCanvas("DifferentSamples4", "DifferentSamples4", 1000, 400); |
1194 | canvas4->Divide(2, 1); |
1195 | |
1196 | TCanvas* canvas3 = new TCanvas("DifferentSamples3", "DifferentSamples3", 1000, 400); |
1197 | canvas3->Divide(2, 1); |
1198 | |
1199 | Float_t min = 1e10; |
1200 | Float_t max = 0; |
1201 | |
1202 | TH1* firstChi = 0; |
1203 | TH1* firstBayesian = 0; |
1204 | |
1205 | TLegend* legend = new TLegend(0.7, 0.7, 1, 1); |
1206 | |
1207 | TFile* file = TFile::Open("DifferentSamples.root", "RECREATE"); |
1208 | file->Close(); |
1209 | |
1210 | for (Int_t i=0; i<8; ++i) |
1211 | { |
1212 | TFile::Open(files[i]); |
1213 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2"); |
1214 | multESD->LoadHistograms("Multiplicity"); |
1215 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
1216 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i)); |
1217 | mc->Sumw2(); |
1218 | |
1219 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
1220 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); |
1221 | mult->DrawComparison(Form("DifferentSamples_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); |
1222 | |
1223 | Float_t chi2MC = 0; |
1224 | Int_t chi2MCLimit = 0; |
1225 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
1226 | fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC); |
1227 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit); |
1228 | min = TMath::Min(min, chi2MC); |
1229 | max = TMath::Max(max, chi2MC); |
1230 | |
1231 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i)); |
1232 | if (!firstChi) |
1233 | firstChi = (TH1*) chi2Result->Clone("firstChi"); |
1234 | |
1235 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100); |
1236 | mult->DrawComparison(Form("DifferentSamples_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
1237 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i)); |
1238 | if (!firstBayesian) |
1239 | firstBayesian = (TH1*) bayesResult->Clone("firstBayesian"); |
1240 | |
1241 | TFile* file = TFile::Open("DifferentSamples.root", "UPDATE"); |
1242 | mc->Write(); |
1243 | chi2Result->Write(); |
1244 | bayesResult->Write(); |
1245 | file->Close(); |
1246 | |
1247 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
1248 | fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC); |
1249 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit); |
1250 | |
1251 | min = TMath::Min(min, chi2MC); |
1252 | max = TMath::Max(max, chi2MC); |
1253 | mc->GetXaxis()->SetRangeUser(0, 150); |
1254 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
1255 | |
1256 | // skip errors for now |
1257 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
1258 | { |
1259 | chi2Result->SetBinError(j, 0); |
1260 | bayesResult->SetBinError(j, 0); |
1261 | } |
1262 | |
1263 | canvas4->cd(1); |
1264 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); |
1265 | tmp->SetTitle("Unfolded/MC;Npart;Ratio"); |
1266 | tmp->Divide(mc); |
1267 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1268 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1269 | tmp->SetLineColor(i+1); |
1270 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1271 | |
1272 | canvas4->cd(2); |
1273 | tmp = (TH1*) bayesResult->Clone("tmp"); |
1274 | tmp->SetTitle("Unfolded/MC;Npart;Ratio"); |
1275 | tmp->Divide(mc); |
1276 | tmp->SetLineColor(i+1); |
1277 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1278 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1279 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1280 | |
1281 | canvas3->cd(1); |
1282 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); |
1283 | tmp->SetTitle("Ratio to first result;Npart;Ratio"); |
1284 | tmp->Divide(firstChi); |
1285 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1286 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1287 | tmp->SetLineColor(i+1); |
1288 | legend->AddEntry(tmp, Form("%d", i)); |
1289 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1290 | |
1291 | canvas3->cd(2); |
1292 | tmp = (TH1*) bayesResult->Clone("tmp"); |
1293 | tmp->SetTitle("Ratio to first result;Npart;Ratio"); |
1294 | tmp->Divide(firstBayesian); |
1295 | tmp->SetLineColor(i+1); |
1296 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1297 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1298 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1299 | |
1300 | if (i < 4) |
1301 | { |
1302 | canvasA->cd(i+1); |
1303 | } |
1304 | else |
1305 | canvasB->cd(i+1-4); |
1306 | |
1307 | mc->SetFillColor(kYellow); |
1308 | mc->DrawCopy(); |
1309 | chi2Result->SetLineColor(kRed); |
1310 | chi2Result->DrawCopy("SAME"); |
1311 | bayesResult->SetLineColor(kBlue); |
1312 | bayesResult->DrawCopy("SAME"); |
1313 | gPad->SetLogy(); |
1314 | |
1315 | if (i < 4) |
1316 | { |
1317 | canvasA->cd(i+5); |
1318 | } |
1319 | else |
1320 | canvasB->cd(i+5-4); |
1321 | |
1322 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); |
1323 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); |
1324 | |
1325 | // skip errors for now |
1326 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
1327 | { |
1328 | chi2Result->SetBinError(j, 0); |
1329 | bayesResult->SetBinError(j, 0); |
1330 | } |
1331 | |
1332 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); |
1333 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
1334 | |
1335 | chi2Result->DrawCopy(""); |
1336 | bayesResult->DrawCopy("SAME"); |
1337 | } |
1338 | |
1339 | canvas3->cd(1); |
1340 | legend->Draw(); |
1341 | |
1342 | canvasA->SaveAs(Form("%s.gif", canvasA->GetName())); |
1343 | canvasB->SaveAs(Form("%s.gif", canvasB->GetName())); |
1344 | |
1345 | TCanvas* canvas2 = new TCanvas("DifferentSamples2", "DifferentSamples2", 800, 400); |
1346 | canvas2->Divide(2, 1); |
1347 | |
1348 | canvas2->cd(1); |
1349 | fitResultsChi2->SetMarkerStyle(20); |
1350 | fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); |
1351 | fitResultsChi2->Draw("AP"); |
1352 | |
1353 | fitResultsBayes->SetMarkerStyle(3); |
1354 | fitResultsBayes->SetMarkerColor(2); |
1355 | fitResultsBayes->Draw("P SAME"); |
1356 | |
1357 | gPad->SetLogy(); |
1358 | |
1359 | canvas2->cd(2); |
1360 | fitResultsChi2Limit->SetMarkerStyle(20); |
1361 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); |
1362 | fitResultsChi2Limit->Draw("AP"); |
1363 | |
1364 | fitResultsBayesLimit->SetMarkerStyle(3); |
1365 | fitResultsBayesLimit->SetMarkerColor(2); |
1366 | fitResultsBayesLimit->Draw("P SAME"); |
1367 | |
1368 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); |
1369 | canvas3->SaveAs(Form("%s.gif", canvas3->GetName())); |
1370 | canvas4->SaveAs(Form("%s.gif", canvas4->GetName())); |
1371 | } |
1372 | |
cfc19dd5 |
1373 | void Merge(Int_t n, const char** files, const char* output) |
1374 | { |
447c325d |
1375 | // const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" }; |
1376 | |
1377 | |
cfc19dd5 |
1378 | gSystem->Load("libPWG0base"); |
1379 | |
1380 | AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n]; |
1381 | TList list; |
1382 | for (Int_t i=0; i<n; ++i) |
1383 | { |
1384 | TString name("Multiplicity"); |
1385 | if (i > 0) |
1386 | name.Form("Multiplicity%d", i); |
1387 | |
1388 | TFile::Open(files[i]); |
1389 | data[i] = new AliMultiplicityCorrection(name, name); |
1390 | data[i]->LoadHistograms("Multiplicity"); |
1391 | if (i > 0) |
1392 | list.Add(data[i]); |
1393 | } |
1394 | |
1395 | data[0]->Merge(&list); |
1396 | |
447c325d |
1397 | //data[0]->DrawHistograms(); |
cfc19dd5 |
1398 | |
1399 | TFile::Open(output, "RECREATE"); |
1400 | data[0]->SaveHistograms(); |
1401 | gFile->Close(); |
1402 | } |
1403 | |
1404 | void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root") |
1405 | { |
1406 | gSystem->Load("libPWG0base"); |
1407 | |
1408 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1409 | |
1410 | TFile::Open(fileName); |
1411 | mult->LoadHistograms("Multiplicity"); |
1412 | |
1413 | TF1* func = 0; |
1414 | |
1415 | if (caseNo >= 4) |
1416 | { |
0b4bfd98 |
1417 | func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 500); |
cfc19dd5 |
1418 | func->SetParNames("scaling", "averagen", "k"); |
1419 | } |
1420 | |
1421 | switch (caseNo) |
1422 | { |
0b4bfd98 |
1423 | case 0: func = new TF1("flat", "1000"); break; |
cfc19dd5 |
1424 | case 1: func = new TF1("flat", "501-x"); break; |
1425 | case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break; |
1426 | case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break; |
dd701109 |
1427 | case 4: func->SetParameters(1e7, 10, 2); break; |
447c325d |
1428 | case 5: func->SetParameters(1, 13, 7); break; |
dd701109 |
1429 | case 6: func->SetParameters(1e7, 30, 4); break; |
0b4bfd98 |
1430 | case 7: func->SetParameters(1e7, 30, 2); break; // *** |
dd701109 |
1431 | case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break; |
cfc19dd5 |
1432 | |
1433 | default: return; |
1434 | } |
1435 | |
447c325d |
1436 | new TCanvas; |
1437 | func->Draw(); |
1438 | |
1439 | mult->SetGenMeasFromFunc(func, 3); |
cfc19dd5 |
1440 | |
dd701109 |
1441 | TFile::Open("out.root", "RECREATE"); |
1442 | mult->SaveHistograms(); |
1443 | |
0b4bfd98 |
1444 | new TCanvas; mult->GetMultiplicityESD(3)->ProjectionY()->DrawCopy(); |
1445 | new TCanvas; mult->GetMultiplicityVtx(3)->ProjectionY()->DrawCopy(); |
1446 | |
dd701109 |
1447 | //mult->ApplyBayesianMethod(2, kFALSE); |
1448 | //mult->ApplyMinuitFit(2, kFALSE); |
cfc19dd5 |
1449 | //mult->ApplyGaussianMethod(2, kFALSE); |
447c325d |
1450 | //mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx); |
dd701109 |
1451 | } |
1452 | |
1453 | void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2) |
1454 | { |
1455 | gSystem->Load("libPWG0base"); |
1456 | |
1457 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1458 | |
1459 | TFile::Open(fileName); |
1460 | mult->LoadHistograms("Multiplicity"); |
1461 | |
1462 | // empty under/overflow bins in x, otherwise Project3D takes them into account |
1463 | TH3* corr = mult->GetCorrelation(corrMatrix); |
1464 | for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j) |
1465 | { |
1466 | for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k) |
1467 | { |
1468 | corr->SetBinContent(0, j, k, 0); |
1469 | corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0); |
1470 | } |
1471 | } |
1472 | |
1473 | TH2* proj = (TH2*) corr->Project3D("zy"); |
1474 | |
1475 | // normalize correction for given nPart |
1476 | for (Int_t i=1; i<=proj->GetNbinsX(); ++i) |
1477 | { |
1478 | Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY()); |
1479 | if (sum <= 0) |
1480 | continue; |
1481 | |
1482 | for (Int_t j=1; j<=proj->GetNbinsY(); ++j) |
1483 | { |
1484 | // npart sum to 1 |
1485 | proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum); |
1486 | proj->SetBinError(i, j, proj->GetBinError(i, j) / sum); |
1487 | } |
1488 | } |
1489 | |
1490 | new TCanvas; |
447c325d |
1491 | proj->DrawCopy("COLZ"); |
dd701109 |
1492 | |
1493 | TH1* scaling = proj->ProjectionY("scaling", 1, 1); |
1494 | scaling->Reset(); |
1495 | scaling->SetMarkerStyle(3); |
1496 | //scaling->GetXaxis()->SetRangeUser(0, 50); |
1497 | TH1* mean = (TH1F*) scaling->Clone("mean"); |
1498 | TH1* width = (TH1F*) scaling->Clone("width"); |
1499 | |
1500 | TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500); |
1501 | lognormal->SetParNames("scaling", "mean", "sigma"); |
1502 | lognormal->SetParameters(1, 1, 1); |
1503 | lognormal->SetParLimits(0, 1, 1); |
1504 | lognormal->SetParLimits(1, 0, 100); |
1505 | lognormal->SetParLimits(2, 1e-3, 1); |
1506 | |
1507 | TF1* nbd = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50); |
1508 | nbd->SetParNames("scaling", "averagen", "k"); |
1509 | nbd->SetParameters(1, 13, 5); |
1510 | nbd->SetParLimits(0, 1, 1); |
1511 | nbd->SetParLimits(1, 1, 100); |
1512 | nbd->SetParLimits(2, 1, 1e8); |
1513 | |
1514 | TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50); |
1515 | poisson->SetParNames("scaling", "k", "deltax"); |
1516 | poisson->SetParameters(1, 1, 0); |
1517 | poisson->SetParLimits(0, 0, 10); |
1518 | poisson->SetParLimits(1, 0.01, 100); |
1519 | poisson->SetParLimits(2, 0, 10); |
1520 | |
1521 | TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50); |
1522 | mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth"); |
1523 | mygaus->SetParameters(1, 0, 1, 1, 0, 1); |
1524 | mygaus->SetParLimits(2, 1e-5, 10); |
1525 | mygaus->SetParLimits(4, 1, 1); |
1526 | mygaus->SetParLimits(5, 1e-5, 10); |
1527 | |
1528 | //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50); |
1529 | TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50); |
1530 | sqrt->SetParNames("ydelta", "exp", "xdelta"); |
1531 | sqrt->SetParameters(0, 0, 1); |
1532 | sqrt->SetParLimits(1, 0, 10); |
1533 | |
1534 | const char* fitWith = "gaus"; |
1535 | |
1536 | for (Int_t i=1; i<=150; ++i) |
1537 | { |
1538 | printf("Fitting %d...\n", i); |
1539 | |
1540 | TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e"); |
447c325d |
1541 | |
dd701109 |
1542 | //hist->GetXaxis()->SetRangeUser(0, 50); |
1543 | //lognormal->SetParameter(0, hist->GetMaximum()); |
1544 | hist->Fit(fitWith, "0 M", ""); |
1545 | |
1546 | TF1* func = hist->GetListOfFunctions()->FindObject(fitWith); |
1547 | |
447c325d |
1548 | if (0 && (i % 5 == 0)) |
dd701109 |
1549 | { |
447c325d |
1550 | pad = new TCanvas; |
dd701109 |
1551 | hist->Draw(); |
1552 | func->Clone()->Draw("SAME"); |
447c325d |
1553 | pad->SetLogy(); |
dd701109 |
1554 | } |
1555 | |
1556 | scaling->Fill(i, func->GetParameter(0)); |
1557 | mean->Fill(i, func->GetParameter(1)); |
1558 | width->Fill(i, func->GetParameter(2)); |
1559 | } |
1560 | |
1561 | TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500); |
1562 | log->SetParameters(0, 1, 1); |
1563 | log->SetParLimits(1, 0, 100); |
1564 | log->SetParLimits(2, 1e-3, 10); |
1565 | |
1566 | TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500); |
1567 | over->SetParameters(0, 1, 0); |
1568 | //over->SetParLimits(0, 0, 100); |
1569 | over->SetParLimits(1, 1e-3, 10); |
1570 | over->SetParLimits(2, 0, 100); |
1571 | |
1572 | c1 = new TCanvas("fitparams", "fitparams", 1200, 400); |
1573 | c1->Divide(3, 1); |
1574 | |
1575 | c1->cd(1); |
1576 | scaling->Draw("P"); |
1577 | |
1578 | //TF1* scalingFit = new TF1("mypol0", "[0]"); |
1579 | TF1* scalingFit = over; |
447c325d |
1580 | scaling->Fit(scalingFit, "", "", 3, 140); |
1581 | scalingFit->SetRange(0, 200); |
1582 | scalingFit->Draw("SAME"); |
dd701109 |
1583 | |
1584 | c1->cd(2); |
1585 | mean->Draw("P"); |
1586 | |
1587 | //TF1* meanFit = log; |
1588 | TF1* meanFit = new TF1("mypol1", "[0]+[1]*x"); |
447c325d |
1589 | mean->Fit(meanFit, "", "", 3, 140); |
1590 | meanFit->SetRange(0, 200); |
1591 | meanFit->Draw("SAME"); |
dd701109 |
1592 | |
1593 | c1->cd(3); |
1594 | width->Draw("P"); |
1595 | |
1596 | //TF1* widthFit = over; |
447c325d |
1597 | TF1* widthFit = new TF1("mypol", "[0]+[1]*TMath::Sqrt([2]*x)"); |
1598 | widthFit->SetParLimits(2, 1e-5, 1e5); |
1599 | width->Fit(widthFit, "", "", 5, 140); |
1600 | widthFit->SetRange(0, 200); |
1601 | widthFit->Draw("SAME"); |
dd701109 |
1602 | |
1603 | // build new correction matrix |
1604 | TH2* new = (TH2*) proj->Clone("new"); |
1605 | new->Reset(); |
1606 | Float_t x, y; |
1607 | for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) |
1608 | { |
1609 | TF1* func = (TF1*) gROOT->FindObject(fitWith); |
1610 | x = new->GetXaxis()->GetBinCenter(i); |
1611 | //if (i == 1) |
1612 | // x = 0.1; |
1613 | x++; |
1614 | func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); |
1615 | printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); |
1616 | |
1617 | for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) |
1618 | { |
447c325d |
1619 | if (i < 11) |
dd701109 |
1620 | { |
1621 | // leave bins 1..20 untouched |
1622 | new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j)); |
1623 | } |
1624 | else |
1625 | { |
1626 | y = new->GetYaxis()->GetBinCenter(j); |
1627 | if (j == 1) |
1628 | y = 0.1; |
1629 | if (func->Eval(y) > 1e-4) |
1630 | new->SetBinContent(i, j, func->Eval(y)); |
1631 | } |
1632 | } |
1633 | } |
1634 | |
1635 | // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0 |
1636 | // we take the values from the old response matrix |
1637 | //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) |
1638 | // new->SetBinContent(i, 1, proj->GetBinContent(i, 1)); |
1639 | |
1640 | //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) |
1641 | // new->SetBinContent(1, j, proj->GetBinContent(1, j)); |
1642 | |
1643 | // normalize correction for given nPart |
1644 | for (Int_t i=1; i<=new->GetNbinsX(); ++i) |
1645 | { |
1646 | Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY()); |
1647 | if (sum <= 0) |
1648 | continue; |
1649 | |
1650 | for (Int_t j=1; j<=new->GetNbinsY(); ++j) |
1651 | { |
1652 | // npart sum to 1 |
1653 | new->SetBinContent(i, j, new->GetBinContent(i, j) / sum); |
1654 | new->SetBinError(i, j, new->GetBinError(i, j) / sum); |
1655 | } |
1656 | } |
1657 | |
1658 | new TCanvas; |
1659 | new->Draw("COLZ"); |
1660 | |
1661 | TH2* diff = (TH2*) new->Clone("diff"); |
1662 | diff->Add(proj, -1); |
1663 | |
1664 | new TCanvas; |
1665 | diff->Draw("COLZ"); |
1666 | diff->SetMinimum(-0.05); |
1667 | diff->SetMaximum(0.05); |
1668 | |
1669 | corr->Reset(); |
1670 | |
1671 | for (Int_t i=1; i<=new->GetNbinsX(); ++i) |
1672 | for (Int_t j=1; j<=new->GetNbinsY(); ++j) |
1673 | corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2, i, j, new->GetBinContent(i, j)); |
1674 | |
1675 | new TCanvas; |
1676 | corr->Project3D("zy")->Draw("COLZ"); |
1677 | |
1678 | TFile::Open("out.root", "RECREATE"); |
1679 | mult->SaveHistograms(); |
1680 | |
1681 | TH1* proj1 = proj->ProjectionY("proj1", 36, 36); |
1682 | TH1* proj2 = new->ProjectionY("proj2", 36, 36); |
1683 | proj2->SetLineColor(2); |
1684 | |
1685 | new TCanvas; |
1686 | proj1->Draw(); |
1687 | proj2->Draw("SAME"); |
cfc19dd5 |
1688 | } |
447c325d |
1689 | |
0b4bfd98 |
1690 | void buildCorrelationMap(const char* fileName = "multiplicityMC_2M.root", Int_t corrMatrix = 3) |
1691 | { |
1692 | gSystem->Load("libPWG0base"); |
1693 | |
1694 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1695 | |
1696 | TFile::Open(fileName); |
1697 | mult->LoadHistograms("Multiplicity"); |
1698 | |
1699 | TH3F* new = mult->GetCorrelation(corrMatrix); |
1700 | new->Reset(); |
1701 | |
1702 | TF1* func = new TF1("func", "gaus(0)"); |
1703 | |
1704 | Int_t vtxBin = new->GetNbinsX() / 2; |
1705 | if (vtxBin == 0) |
1706 | vtxBin = 1; |
1707 | |
1708 | Float_t sigma = 2; |
1709 | for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1) |
1710 | { |
1711 | Float_t x = new->GetYaxis()->GetBinCenter(i); |
1712 | func->SetParameters(1, x * 0.8, sigma); |
1713 | //func->SetParameters(1, x, sigma); |
1714 | |
1715 | for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1) |
1716 | { |
1717 | Float_t y = new->GetYaxis()->GetBinCenter(j); |
1718 | |
1719 | // cut at 1 sigma |
1720 | if (TMath::Abs(y-x*0.8) < sigma) |
1721 | new->SetBinContent(vtxBin, i, j, func->Eval(y)); |
1722 | |
1723 | // test only bin 40 has smearing |
1724 | //if (x != 40) |
1725 | // new->SetBinContent(vtxBin, i, j, (i == j)); |
1726 | } |
1727 | } |
1728 | |
1729 | new TCanvas; |
1730 | new->Project3D("zy")->DrawCopy("COLZ"); |
1731 | |
1732 | TFile* file = TFile::Open("out.root", "RECREATE"); |
1733 | mult->SetCorrelation(corrMatrix, new); |
1734 | mult->SaveHistograms(); |
1735 | file->Close(); |
1736 | } |
1737 | |
447c325d |
1738 | void GetCrossSections(const char* fileName) |
1739 | { |
1740 | gSystem->Load("libPWG0base"); |
1741 | |
1742 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1743 | |
1744 | TFile::Open(fileName); |
1745 | mult->LoadHistograms("Multiplicity"); |
1746 | |
1747 | TH1* xSection2 = mult->GetCorrelation(3)->Project3D("y")->Clone("xSection2"); |
1748 | xSection2->Sumw2(); |
1749 | xSection2->Scale(1.0 / xSection2->Integral()); |
1750 | |
1751 | TH1* xSection15 = mult->GetCorrelation(2)->Project3D("y")->Clone("xSection15"); |
1752 | xSection15->Sumw2(); |
1753 | xSection15->Scale(1.0 / xSection15->Integral()); |
1754 | |
1755 | TFile::Open("crosssection.root", "RECREATE"); |
1756 | xSection2->Write(); |
1757 | xSection15->Write(); |
1758 | gFile->Close(); |
1759 | } |
0b4bfd98 |
1760 | |
44466df2 |
1761 | void AnalyzeSpeciesTree(const char* fileName) |
1762 | { |
1763 | // |
1764 | // prints statistics about fParticleSpecies |
1765 | // |
1766 | |
1767 | gSystem->Load("libPWG0base"); |
1768 | |
1769 | TFile::Open(fileName); |
1770 | TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies"); |
1771 | |
1772 | const Int_t nFields = 8; |
1773 | Long_t totals[8]; |
1774 | for (Int_t i=0; i<nFields; i++) |
1775 | totals[i] = 0; |
1776 | |
1777 | for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++) |
1778 | { |
1779 | fParticleSpecies->GetEvent(i); |
1780 | |
1781 | Float_t* f = fParticleSpecies->GetArgs(); |
1782 | |
1783 | for (Int_t j=0; j<nFields; j++) |
1784 | totals[j] += f[j+1]; |
1785 | } |
1786 | |
1787 | for (Int_t i=0; i<nFields; i++) |
1788 | Printf("%d --> %ld", i, totals[i]); |
1789 | } |
1790 | |
0b4bfd98 |
1791 | void BuildResponseFromTree(const char* fileName, const char* target) |
1792 | { |
1793 | // |
1794 | // builds several response matrices with different particle ratios (systematic study) |
1795 | // |
1796 | |
1797 | gSystem->Load("libPWG0base"); |
1798 | |
1799 | TFile::Open(fileName); |
1800 | TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies"); |
1801 | |
1802 | TFile* file = TFile::Open(target, "RECREATE"); |
1803 | file->Close(); |
1804 | |
1805 | Int_t tracks = 0; // control variables |
1806 | Int_t noLabel = 0; |
6d81c2de |
1807 | Int_t secondaries = 0; |
0b4bfd98 |
1808 | Int_t doubleCount = 0; |
1809 | |
44466df2 |
1810 | for (Int_t num = 0; num < 7; num++) |
0b4bfd98 |
1811 | { |
1812 | AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(Form("Multiplicity_%d", num), Form("Multiplicity_%d", num)); |
1813 | |
1814 | Float_t ratio[4]; // pi, K, p, other |
1815 | for (Int_t i = 0; i < 4; i++) |
1816 | ratio[i] = 1; |
1817 | |
1818 | switch (num) |
1819 | { |
1820 | case 1 : ratio[1] = 0.5; break; |
1821 | case 2 : ratio[2] = 0.5; break; |
1822 | case 3 : ratio[1] = 1.5; break; |
1823 | case 4 : ratio[2] = 1.5; break; |
1824 | case 5 : ratio[1] = 0.5; ratio[2] = 0.5; break; |
1825 | case 6 : ratio[1] = 1.5; ratio[2] = 1.5; break; |
1826 | } |
1827 | |
1828 | for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++) |
1829 | { |
1830 | fParticleSpecies->GetEvent(i); |
1831 | |
1832 | Float_t* f = fParticleSpecies->GetArgs(); |
1833 | |
1834 | Float_t gene = 0; |
1835 | Float_t meas = 0; |
1836 | |
1837 | for (Int_t j = 0; j < 4; j++) |
1838 | { |
1839 | gene += ratio[j] * f[j+1]; |
1840 | meas += ratio[j] * f[j+1+4]; |
1841 | tracks += f[j+1+4]; |
1842 | } |
1843 | |
1844 | // add the ones w/o label |
1845 | tracks += f[9]; |
1846 | noLabel += f[9]; |
1847 | |
6d81c2de |
1848 | // secondaries are already part of meas! |
1849 | secondaries += f[10]; |
1850 | |
1851 | // double counted are already part of meas! |
1852 | doubleCount += f[11]; |
0b4bfd98 |
1853 | |
6d81c2de |
1854 | // ones w/o label are added without weighting to allow comparison to default analysis. however this is only valid when their fraction is low enough! |
0b4bfd98 |
1855 | meas += f[9]; |
0b4bfd98 |
1856 | |
1857 | //Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]); |
1858 | |
6d81c2de |
1859 | fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, 0, meas, meas, meas, meas); |
1860 | fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, gene, gene, gene, gene, 0); |
1861 | fMultiplicity->FillMeasured(f[0], meas, meas, meas, meas); |
0b4bfd98 |
1862 | } |
1863 | |
1864 | //fMultiplicity->DrawHistograms(); |
1865 | |
1866 | TFile* file = TFile::Open(target, "UPDATE"); |
1867 | fMultiplicity->SaveHistograms(); |
1868 | file->Close(); |
1869 | |
1870 | if (num == 0) |
6d81c2de |
1871 | { |
1872 | Printf("%d total tracks, %d w/o label = %.2f %%, %d double counted = %.2f %%, secondaries = %.2f %%", tracks, noLabel, 100.0 * noLabel / tracks, doubleCount, 100.0 * doubleCount / tracks, 100.0 * secondaries / tracks); |
1873 | if ((Float_t) noLabel / tracks > 0.02) |
1874 | Printf("WARNING: More than 2%% of tracks without label, this might bias the study!"); |
1875 | } |
0b4bfd98 |
1876 | } |
1877 | } |
1878 | |
1879 | void MergeModifyCrossSection(const char* output) |
1880 | { |
1881 | const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root" }; |
1882 | |
1883 | gSystem->Load("libPWG0base"); |
1884 | |
1885 | TFile::Open(output, "RECREATE"); |
1886 | gFile->Close(); |
1887 | |
1888 | for (Int_t num=0; num<7; ++num) |
1889 | { |
1890 | AliMultiplicityCorrection* data[3]; |
1891 | TList list; |
1892 | |
1893 | Float_t ratio[3]; |
1894 | switch (num) |
1895 | { |
1896 | case 0: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.0; break; |
1897 | case 1: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.0; break; |
1898 | case 2: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 1.0; break; |
1899 | case 3: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.5; break; |
1900 | case 4: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 0.5; break; |
1901 | case 5: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.5; break; |
1902 | case 6: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 0.5; break; |
1903 | default: return; |
1904 | } |
1905 | |
1906 | for (Int_t i=0; i<3; ++i) |
1907 | { |
1908 | TString name; |
1909 | name.Form("Multiplicity_%d", num); |
1910 | if (i > 0) |
1911 | name.Form("Multiplicity_%d_%d", num, i); |
1912 | |
1913 | TFile::Open(files[i]); |
1914 | data[i] = new AliMultiplicityCorrection(name, name); |
1915 | data[i]->LoadHistograms("Multiplicity"); |
1916 | |
1917 | // modify x-section |
1918 | for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++) |
1919 | { |
1920 | data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]); |
1921 | data[i]->GetMultiplicityMB(j)->Scale(ratio[i]); |
1922 | data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]); |
1923 | } |
1924 | |
1925 | for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++) |
1926 | data[i]->GetMultiplicityESD(j)->Scale(ratio[i]); |
1927 | |
1928 | for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++) |
1929 | data[i]->GetCorrelation(j)->Scale(ratio[i]); |
1930 | |
1931 | if (i > 0) |
1932 | list.Add(data[i]); |
1933 | } |
1934 | |
1935 | printf("Case %d, %s: Entries in response matrix 3: ND: %.2f SD: %.2f DD: %.2f", num, data[0]->GetName(), data[0]->GetCorrelation(3)->Integral(), data[1]->GetCorrelation(3)->Integral(), data[2]->GetCorrelation(3)->Integral()); |
1936 | |
1937 | data[0]->Merge(&list); |
1938 | |
1939 | Printf(" Total: %.2f", data[0]->GetCorrelation(3)->Integral()); |
1940 | |
1941 | TFile::Open(output, "UPDATE"); |
1942 | data[0]->SaveHistograms(); |
1943 | gFile->Close(); |
1944 | |
1945 | list.Clear(); |
1946 | |
1947 | for (Int_t i=0; i<3; ++i) |
1948 | delete data[i]; |
1949 | } |
1950 | } |
44466df2 |
1951 | |
1952 | void Rebin(const char* fileName = "multiplicityMC_3M.root", Int_t corrMatrix = 3) |
1953 | { |
1954 | // rebins MC axis of correlation map, MC and histogram for corrected (for evaluation of effect of regularization) |
1955 | // rebin does not exist for 3D hists, so we convert to 2D and then back to 3D (loosing the vertex information) |
1956 | |
1957 | Printf("WARNING: Vertex information is lost in this process. Use result only for evaluation of errors."); |
1958 | |
1959 | gSystem->Load("libPWG0base"); |
1960 | |
1961 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1962 | |
1963 | TFile::Open(fileName); |
1964 | mult->LoadHistograms("Multiplicity"); |
1965 | |
1966 | // rebin correlation |
1967 | TH3* old = mult->GetCorrelation(corrMatrix); |
1968 | |
1969 | // empty under/overflow bins in x, otherwise Project3D takes them into account |
1970 | for (Int_t y=1; y<=old->GetYaxis()->GetNbins(); ++y) |
1971 | { |
1972 | for (Int_t z=1; z<=old->GetZaxis()->GetNbins(); ++z) |
1973 | { |
1974 | old->SetBinContent(0, y, z, 0); |
1975 | old->SetBinContent(old->GetXaxis()->GetNbins()+1, y, z, 0); |
1976 | } |
1977 | } |
1978 | |
1979 | TH2* response = (TH2*) old->Project3D("zy"); |
1980 | response->RebinX(2); |
1981 | |
1982 | TH3F* new = new TH3F(old->GetName(), old->GetTitle(), |
1983 | old->GetXaxis()->GetNbins(), old->GetXaxis()->GetBinLowEdge(1), old->GetXaxis()->GetBinUpEdge(old->GetXaxis()->GetNbins()), |
1984 | old->GetYaxis()->GetNbins() / 2, old->GetYaxis()->GetBinLowEdge(1), old->GetYaxis()->GetBinUpEdge(old->GetYaxis()->GetNbins()), |
1985 | old->GetZaxis()->GetNbins(), old->GetZaxis()->GetBinLowEdge(1), old->GetZaxis()->GetBinUpEdge(old->GetZaxis()->GetNbins())); |
1986 | new->Reset(); |
1987 | |
1988 | Int_t vtxBin = new->GetNbinsX() / 2; |
1989 | if (vtxBin == 0) |
1990 | vtxBin = 1; |
1991 | |
1992 | for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1) |
1993 | for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1) |
1994 | new->SetBinContent(vtxBin, i, j, response->GetBinContent(i, j)); |
1995 | |
1996 | // rebin MC + hist for corrected |
1997 | for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kTrVtx; eventType <= AliMultiplicityCorrection::kINEL; eventType++) |
1998 | mult->GetMultiplicityMC(corrMatrix, eventType)->RebinY(2); |
1999 | |
2000 | mult->GetMultiplicityESDCorrected(corrMatrix)->Rebin(2); |
2001 | |
2002 | // recreate measured from correlation matrix to get rid of vertex shift effect |
2003 | TH2* newMeasured = (TH2*) old->Project3D("zx"); |
2004 | TH2* esd = mult->GetMultiplicityESD(corrMatrix); |
2005 | esd->Reset(); |
2006 | |
2007 | // transfer from TH2D to TH2F |
2008 | for (Int_t i=0; i<=new->GetXaxis()->GetNbins()+1; i+=1) |
2009 | for (Int_t j=0; j<=new->GetYaxis()->GetNbins()+1; j+=1) |
2010 | esd->SetBinContent(i, j, newMeasured->GetBinContent(i, j)); |
2011 | |
2012 | new TCanvas; |
2013 | new->Project3D("zy")->DrawCopy("COLZ"); |
2014 | |
2015 | TFile* file = TFile::Open("out.root", "RECREATE"); |
2016 | mult->SetCorrelation(corrMatrix, new); |
2017 | mult->SaveHistograms(); |
2018 | file->Close(); |
2019 | } |
2020 | |
2021 | void EvaluateRegularizationEffect(Int_t step, const char* fileNameRebinned = "multiplicityMC_3M_rebinned.root", const char* fileNameNormal = "multiplicityMC_3M.root", Int_t histID = 3) |
2022 | { |
2023 | // due to some static members in AliMultiplicityCorrection, the session has to be restarted after changing the number of parameters, to be fixed |
2024 | // that is why this is done in 2 steps |
2025 | |
2026 | gSystem->Load("libPWG0base"); |
2027 | |
2028 | Bool_t fullPhaseSpace = kFALSE; |
2029 | |
2030 | if (step == 1) |
2031 | { |
2032 | // first step: unfold without regularization and rebinned histogram ("N=M") |
2033 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
2034 | TFile::Open(fileNameRebinned); |
2035 | mult->LoadHistograms(); |
2036 | |
2037 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); |
2038 | mult->SetCreateBigBin(kFALSE); |
2039 | |
2040 | mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); |
2041 | mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, mult->GetMultiplicityVtx(histID)->ProjectionY("mymchist")); |
2042 | |
2043 | TFile* file = TFile::Open("EvaluateRegularizationEffect1.root", "RECREATE"); |
2044 | mult->SaveHistograms(); |
2045 | file->Close(); |
2046 | } |
2047 | else if (step == 2) |
2048 | { |
2049 | // second step: unfold with regularization and normal histogram |
2050 | AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
2051 | TFile::Open(fileNameNormal); |
2052 | mult2->LoadHistograms(); |
2053 | |
2054 | mult2->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4); |
2055 | mult2->SetCreateBigBin(kTRUE); |
2056 | mult2->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); |
2057 | mult2->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, mult2->GetMultiplicityVtx(histID)->ProjectionY("mymchist")); |
2058 | |
2059 | TH1* result2 = mult2->GetMultiplicityESDCorrected(histID); |
2060 | |
2061 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
2062 | TFile* file = TFile::Open("EvaluateRegularizationEffect1.root"); |
2063 | mult->LoadHistograms(); |
2064 | |
2065 | TH1* result1 = mult->GetMultiplicityESDCorrected(histID); |
2066 | |
2067 | // compare results |
2068 | TCanvas* canvas = new TCanvas("EvaluateRegularizationEffect", "EvaluateRegularizationEffect", 1000, 800); |
2069 | canvas->Divide(2, 2); |
2070 | |
2071 | canvas->cd(1); |
2072 | result1->SetLineColor(1); |
2073 | result1->DrawCopy(); |
2074 | result2->SetLineColor(2); |
2075 | result2->DrawCopy("SAME"); |
2076 | gPad->SetLogy(); |
2077 | |
2078 | result2->Rebin(2); |
2079 | result1->Scale(1.0 / result1->Integral()); |
2080 | result2->Scale(1.0 / result2->Integral()); |
2081 | |
2082 | canvas->cd(2); |
2083 | result1->DrawCopy(); |
2084 | result2->DrawCopy("SAME"); |
2085 | gPad->SetLogy(); |
2086 | |
2087 | TH1* diff = (TH1*) result1->Clone("diff"); |
2088 | diff->Add(result2, -1); |
2089 | |
2090 | canvas->cd(3); |
2091 | diff->DrawCopy("HIST"); |
2092 | |
2093 | canvas->cd(4); |
2094 | diff->Divide(result1); |
2095 | diff->GetYaxis()->SetRangeUser(-0.3, 0.3); |
2096 | diff->DrawCopy("HIST"); |
2097 | |
2098 | Double_t chi2 = 0; |
2099 | for (Int_t i=1; i<=diff->GetNbinsX(); i++) |
2100 | chi2 += diff->GetBinContent(i) * diff->GetBinContent(i); |
2101 | |
2102 | Printf("Chi2 is %e", chi2); |
2103 | |
2104 | canvas->SaveAs(Form("%s.eps", canvas->GetName())); |
2105 | } |
2106 | } |