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