0ab29cfa |
1 | /* $Id$ */ |
2 | |
3 | // |
4 | // script to run the AliMultiplicityESDSelector |
5 | // |
6 | |
7 | #include "../CreateESDChain.C" |
0bd1f8a0 |
8 | #include "../PWG0Helper.C" |
0ab29cfa |
9 | |
03244034 |
10 | void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = "") |
0ab29cfa |
11 | { |
0bd1f8a0 |
12 | if (aProof) |
cfc19dd5 |
13 | connectProof("jgrosseo@lxb6046"); |
0bd1f8a0 |
14 | |
15 | TString libraries("libEG;libGeom;libESD;libPWG0base"); |
16 | TString packages("PWG0base"); |
10ebe68d |
17 | |
0ab29cfa |
18 | if (aMC != kFALSE) |
0bd1f8a0 |
19 | { |
20 | libraries += ";libVMC;libMinuit;libSTEER;libPWG0dep;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6"; |
21 | packages += ";PWG0dep"; |
22 | } |
23 | |
24 | if (!prepareQuery(libraries, packages, kTRUE)) |
25 | return; |
0ab29cfa |
26 | |
0ab29cfa |
27 | gROOT->ProcessLine(".L CreateCuts.C"); |
28 | gROOT->ProcessLine(".L drawPlots.C"); |
29 | |
0ab29cfa |
30 | // selection of esd tracks |
31 | AliESDtrackCuts* esdTrackCuts = CreateTrackCuts(); |
32 | if (!esdTrackCuts) |
33 | { |
34 | printf("ERROR: esdTrackCuts could not be created\n"); |
35 | return; |
36 | } |
37 | |
0bd1f8a0 |
38 | TList inputList; |
39 | inputList.Add(esdTrackCuts); |
40 | |
cfc19dd5 |
41 | TChain* chain = CreateESDChain(data, nRuns, offset, kFALSE, kFALSE); |
0ab29cfa |
42 | |
43 | TString selectorName = ((aMC == kFALSE) ? "AliMultiplicityESDSelector" : "AliMultiplicityMCSelector"); |
44 | AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo); |
45 | |
0bd1f8a0 |
46 | selectorName += ".cxx+"; |
0ab29cfa |
47 | |
48 | if (aDebug != kFALSE) |
49 | selectorName += "g"; |
50 | |
03244034 |
51 | Int_t result = executeQuery(chain, &inputList, selectorName, option); |
0ab29cfa |
52 | |
53 | if (result != 0) |
54 | { |
55 | printf("ERROR: Executing process failed with %d.\n", result); |
56 | return; |
57 | } |
0a173978 |
58 | } |
0ab29cfa |
59 | |
0a173978 |
60 | void draw(const char* fileName = "multiplicityMC.root") |
61 | { |
62 | gSystem->Load("libPWG0base"); |
63 | |
64 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
65 | |
66 | TFile::Open(fileName); |
67 | mult->LoadHistograms("Multiplicity"); |
68 | |
69 | mult->DrawHistograms(); |
0ab29cfa |
70 | } |
71 | |
cfc19dd5 |
72 | void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0) |
0a173978 |
73 | { |
74 | gSystem->Load("libPWG0base"); |
75 | |
76 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
77 | |
78 | TFile::Open(fileName); |
79 | mult->LoadHistograms("Multiplicity"); |
80 | |
dd701109 |
81 | //mult->ApplyLaszloMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx); |
82 | |
83 | //return; |
84 | |
85 | |
86 | mult->ApplyBayesianMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx); |
87 | mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY()); |
88 | |
89 | return; |
90 | |
91 | TStopwatch timer; |
cfc19dd5 |
92 | |
dd701109 |
93 | timer.Start(); |
94 | |
95 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1); |
96 | mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx); |
97 | mult->DrawComparison("MinuitChi2", hist, kFALSE, kFALSE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY()); |
98 | |
99 | timer.Stop(); |
100 | timer.Print(); |
101 | |
102 | return 0; |
cfc19dd5 |
103 | |
104 | //mult->ApplyGaussianMethod(hist, kFALSE); |
9ca6be09 |
105 | |
dd701109 |
106 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); |
107 | mult->ApplyNBDFit(hist, kFALSE); |
108 | mult->DrawComparison("NBDChi2Fit", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY()); |
109 | |
9ca6be09 |
110 | return mult; |
111 | } |
112 | |
dd701109 |
113 | void* fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Int_t histID = 2) |
9ca6be09 |
114 | { |
115 | gSystem->Load("libPWG0base"); |
116 | |
117 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
118 | |
119 | TFile::Open(fileNameMC); |
120 | mult->LoadHistograms("Multiplicity"); |
121 | |
122 | TFile::Open(fileNameESD); |
cfc19dd5 |
123 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); |
124 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID)); |
dd701109 |
125 | //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityMB%d", histID)); |
9ca6be09 |
126 | |
cfc19dd5 |
127 | mult->SetMultiplicityESD(histID, hist); |
9ca6be09 |
128 | |
dd701109 |
129 | mult->SetMultiplicityVtx(histID, hist2); |
130 | mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
131 | return; |
132 | |
133 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1.1); |
134 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
135 | mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY()); |
136 | |
137 | return; |
138 | |
cfc19dd5 |
139 | //mult->ApplyGaussianMethod(histID, kFALSE); |
cfc19dd5 |
140 | |
dd701109 |
141 | for (Float_t f=0.1; f<=0.11; f+=0.05) |
142 | { |
143 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, f); |
144 | mult->DrawComparison(Form("Bayesian_%f", f), histID, kFALSE, kTRUE, hist2->ProjectionY()); |
145 | } |
146 | |
147 | //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e7); |
148 | //mult->ApplyMinuitFit(histID, kFALSE); |
149 | //mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY()); |
150 | |
0a173978 |
151 | |
152 | return mult; |
153 | } |
cfc19dd5 |
154 | |
155 | const char* GetRegName(Int_t type) |
156 | { |
157 | switch (type) |
158 | { |
159 | case AliMultiplicityCorrection::kNone: return "None"; break; |
160 | case AliMultiplicityCorrection::kPol0: return "Pol0"; break; |
161 | case AliMultiplicityCorrection::kPol1: return "Pol1"; break; |
162 | case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break; |
163 | case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break; |
dd701109 |
164 | case AliMultiplicityCorrection::kTest : return "Test"; break; |
165 | } |
166 | return 0; |
167 | } |
168 | |
169 | const char* GetEventTypeName(Int_t type) |
170 | { |
171 | switch (type) |
172 | { |
173 | case AliMultiplicityCorrection::kTrVtx: return "trigger, vertex"; break; |
174 | case AliMultiplicityCorrection::kMB: return "minimum bias"; break; |
175 | case AliMultiplicityCorrection::kINEL: return "inelastic"; break; |
cfc19dd5 |
176 | } |
177 | return 0; |
178 | } |
179 | |
dd701109 |
180 | void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) |
cfc19dd5 |
181 | { |
dd701109 |
182 | gSystem->mkdir(targetDir); |
183 | |
184 | gSystem->Load("libPWG0base"); |
185 | |
186 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
187 | TFile::Open(fileNameMC); |
188 | mult->LoadHistograms("Multiplicity"); |
189 | |
190 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
191 | TFile::Open(fileNameESD); |
192 | multESD->LoadHistograms("Multiplicity"); |
193 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
194 | |
195 | TCanvas* canvas = new TCanvas("EvaluateBayesianMethod", "EvaluateBayesianMethod", 800, 600); |
196 | TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); |
197 | legend->SetFillColor(0); |
198 | |
199 | Float_t min = 1e10; |
200 | Float_t max = 0; |
201 | |
202 | TGraph* first = 0; |
203 | Int_t count = 0; // just to order the saved images... |
204 | |
205 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) |
206 | { |
207 | TGraph* fitResultsMC = new TGraph; |
208 | fitResultsMC->SetTitle(";Weight Parameter"); |
209 | TGraph* fitResultsRes = new TGraph; |
210 | fitResultsRes->SetTitle(";Weight Parameter"); |
211 | |
212 | fitResultsMC->SetFillColor(0); |
213 | fitResultsRes->SetFillColor(0); |
214 | fitResultsMC->SetMarkerStyle(20+type); |
215 | fitResultsRes->SetMarkerStyle(24+type); |
216 | fitResultsRes->SetMarkerColor(kRed); |
217 | fitResultsRes->SetLineColor(kRed); |
218 | |
219 | legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type))); |
220 | legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type))); |
221 | |
222 | for (Float_t weight = 0; weight < 0.301; weight += 0.02) |
223 | { |
224 | Float_t chi2MC = 0; |
225 | Float_t residuals = 0; |
226 | |
227 | mult->ApplyBayesianMethod(histID, kFALSE, type, weight); |
228 | mult->DrawComparison(Form("%s/Bayesian_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); |
229 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
230 | |
231 | fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC); |
232 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); |
233 | |
234 | min = TMath::Min(min, TMath::Min(chi2MC, residuals)); |
235 | max = TMath::Max(max, TMath::Max(chi2MC, residuals)); |
236 | } |
237 | |
238 | fitResultsMC->Print(); |
239 | fitResultsRes->Print(); |
240 | |
241 | canvas->cd(); |
242 | fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME")); |
243 | fitResultsRes->Draw("SAME CP"); |
244 | |
245 | if (first == 0) |
246 | first = fitResultsMC; |
247 | } |
248 | |
249 | gPad->SetLogy(); |
250 | printf("min = %f, max = %f\n", min, max); |
251 | if (min <= 0) |
252 | min = 1e-5; |
253 | first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5); |
254 | |
255 | legend->Draw(); |
256 | |
257 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
258 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); |
259 | } |
260 | |
261 | void EvaluateBayesianMethodIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) |
262 | { |
263 | gSystem->mkdir(targetDir); |
264 | |
265 | gSystem->Load("libPWG0base"); |
266 | |
267 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
268 | TFile::Open(fileNameMC); |
269 | mult->LoadHistograms("Multiplicity"); |
270 | |
271 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
272 | TFile::Open(fileNameESD); |
273 | multESD->LoadHistograms("Multiplicity"); |
274 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
275 | |
276 | TCanvas* canvas = new TCanvas("EvaluateBayesianMethodIterations", "EvaluateBayesianMethodIterations", 800, 600); |
277 | TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); |
278 | legend->SetFillColor(0); |
279 | |
280 | Float_t min = 1e10; |
281 | Float_t max = 0; |
282 | |
283 | TGraph* first = 0; |
284 | Int_t count = 0; // just to order the saved images... |
285 | |
286 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) |
287 | { |
288 | TGraph* fitResultsMC = new TGraph; |
289 | fitResultsMC->SetTitle(";Iterations"); |
290 | TGraph* fitResultsRes = new TGraph; |
291 | fitResultsRes->SetTitle(";Iterations"); |
292 | |
293 | fitResultsMC->SetFillColor(0); |
294 | fitResultsRes->SetFillColor(0); |
295 | fitResultsMC->SetMarkerStyle(20+type); |
296 | fitResultsRes->SetMarkerStyle(24+type); |
297 | fitResultsRes->SetMarkerColor(kRed); |
298 | fitResultsRes->SetLineColor(kRed); |
299 | |
300 | legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type))); |
301 | legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type))); |
302 | |
303 | for (Int_t iter = 5; iter <= 50; iter += 5) |
304 | { |
305 | Float_t chi2MC = 0; |
306 | Float_t residuals = 0; |
307 | |
308 | mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1, iter); |
309 | mult->DrawComparison(Form("%s/BayesianIter_%02d_%d_%d", targetDir, count++, type, iter), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); |
310 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
311 | |
312 | fitResultsMC->SetPoint(fitResultsMC->GetN(), iter, chi2MC); |
313 | fitResultsRes->SetPoint(fitResultsRes->GetN(), iter, residuals); |
314 | |
315 | min = TMath::Min(min, TMath::Min(chi2MC, residuals)); |
316 | max = TMath::Max(max, TMath::Max(chi2MC, residuals)); |
317 | } |
318 | |
319 | fitResultsMC->Print(); |
320 | fitResultsRes->Print(); |
321 | |
322 | canvas->cd(); |
323 | fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME")); |
324 | fitResultsRes->Draw("SAME CP"); |
325 | |
326 | if (first == 0) |
327 | first = fitResultsMC; |
328 | } |
329 | |
330 | gPad->SetLogy(); |
331 | printf("min = %f, max = %f\n", min, max); |
332 | if (min <= 0) |
333 | min = 1e-5; |
334 | first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5); |
335 | |
336 | legend->Draw(); |
337 | |
338 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
339 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); |
340 | } |
341 | |
342 | void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) |
343 | { |
344 | gSystem->mkdir(targetDir); |
345 | |
cfc19dd5 |
346 | gSystem->Load("libPWG0base"); |
347 | |
348 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
349 | |
350 | TFile::Open(fileNameMC); |
351 | mult->LoadHistograms("Multiplicity"); |
352 | |
353 | TFile::Open(fileNameESD); |
354 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); |
355 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID)); |
356 | |
357 | mult->SetMultiplicityESD(histID, hist); |
358 | |
359 | TCanvas* canvas = new TCanvas("EvaluateChi2Method", "EvaluateChi2Method", 800, 600); |
dd701109 |
360 | TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); |
cfc19dd5 |
361 | legend->SetFillColor(0); |
362 | |
363 | Float_t min = 1e10; |
364 | Float_t max = 0; |
365 | |
366 | TGraph* first = 0; |
dd701109 |
367 | Int_t count = 0; // just to order the saved images... |
368 | |
369 | Bool_t firstLoop = kTRUE; |
cfc19dd5 |
370 | |
371 | for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type) |
dd701109 |
372 | //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type) |
cfc19dd5 |
373 | { |
374 | TGraph* fitResultsMC = new TGraph; |
375 | fitResultsMC->SetTitle(";Weight Parameter"); |
376 | TGraph* fitResultsRes = new TGraph; |
377 | fitResultsRes->SetTitle(";Weight Parameter"); |
378 | |
379 | fitResultsMC->SetFillColor(0); |
380 | fitResultsRes->SetFillColor(0); |
381 | fitResultsMC->SetMarkerStyle(19+type); |
dd701109 |
382 | fitResultsRes->SetMarkerStyle(23+type); |
cfc19dd5 |
383 | fitResultsRes->SetMarkerColor(kRed); |
384 | fitResultsRes->SetLineColor(kRed); |
385 | |
386 | legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetRegName(type))); |
387 | legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetRegName(type))); |
388 | |
389 | if (first == 0) |
390 | first = fitResultsMC; |
391 | |
dd701109 |
392 | for (Float_t weight = 1e-4; weight < 2e4; weight *= 10) |
393 | //for (Float_t weight = 0.1; weight < 10; weight *= TMath::Sqrt(TMath::Sqrt(10))) |
cfc19dd5 |
394 | { |
395 | Float_t chi2MC = 0; |
396 | Float_t residuals = 0; |
397 | |
398 | mult->SetRegularizationParameters(type, weight); |
dd701109 |
399 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
400 | mult->DrawComparison(Form("%s/MinuitChi2_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY()); |
401 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
cfc19dd5 |
402 | |
403 | fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC); |
404 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); |
405 | |
406 | min = TMath::Min(min, TMath::Min(chi2MC, residuals)); |
407 | max = TMath::Max(max, TMath::Max(chi2MC, residuals)); |
408 | } |
409 | |
410 | fitResultsMC->Print(); |
411 | fitResultsRes->Print(); |
412 | |
413 | canvas->cd(); |
dd701109 |
414 | fitResultsMC->Draw(Form("%s CP", (firstLoop) ? "A" : "SAME")); |
cfc19dd5 |
415 | fitResultsRes->Draw("SAME CP"); |
dd701109 |
416 | |
417 | firstLoop = kFALSE; |
cfc19dd5 |
418 | } |
419 | |
420 | gPad->SetLogx(); |
421 | gPad->SetLogy(); |
422 | printf("min = %f, max = %f\n", min, max); |
423 | if (min <= 0) |
424 | min = 1e-5; |
425 | first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5); |
426 | |
427 | legend->Draw(); |
428 | |
dd701109 |
429 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
430 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); |
431 | } |
432 | |
433 | void EvaluateChi2MethodAll() |
434 | { |
435 | EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); |
436 | EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); |
437 | EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); |
438 | EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); |
439 | EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); |
440 | } |
441 | |
442 | void EvaluateBayesianMethodAll() |
443 | { |
444 | EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); |
445 | EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); |
446 | EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); |
447 | EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); |
448 | EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); |
449 | } |
450 | |
451 | void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) |
452 | { |
453 | gSystem->mkdir(targetDir); |
454 | |
455 | gSystem->Load("libPWG0base"); |
456 | |
457 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
458 | |
459 | TFile::Open(fileNameMC); |
460 | mult->LoadHistograms("Multiplicity"); |
461 | |
462 | TFile::Open(fileNameESD); |
463 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
464 | multESD->LoadHistograms("Multiplicity"); |
465 | |
466 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
467 | |
468 | TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 800); |
469 | canvas->Divide(3, 2); |
470 | |
471 | Int_t count = 0; |
472 | |
473 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) |
474 | { |
475 | TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY(); |
476 | |
477 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3); |
478 | mult->ApplyMinuitFit(histID, kFALSE, type); |
479 | mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc); |
480 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result"); |
481 | |
482 | mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1); |
483 | mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc); |
484 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult"); |
485 | |
486 | mc->GetXaxis()->SetRangeUser(0, 150); |
487 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
488 | |
489 | // skip errors for now |
490 | for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i) |
491 | { |
492 | chi2Result->SetBinError(i, 0); |
493 | bayesResult->SetBinError(i, 0); |
494 | } |
495 | |
496 | canvas->cd(++count); |
497 | mc->SetFillColor(kYellow); |
498 | mc->DrawCopy(); |
499 | chi2Result->SetLineColor(kRed); |
500 | chi2Result->DrawCopy("SAME"); |
501 | bayesResult->SetLineColor(kBlue); |
502 | bayesResult->DrawCopy("SAME"); |
503 | gPad->SetLogy(); |
504 | |
505 | canvas->cd(count + 3); |
506 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); |
507 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); |
508 | |
509 | // skip errors for now |
510 | for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i) |
511 | { |
512 | chi2Result->SetBinError(i, 0); |
513 | bayesResult->SetBinError(i, 0); |
514 | } |
515 | |
516 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
517 | |
518 | chi2Result->DrawCopy(""); |
519 | bayesResult->DrawCopy("SAME"); |
520 | } |
521 | |
522 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
523 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); |
524 | } |
525 | |
526 | void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2) |
527 | { |
528 | gSystem->Load("libPWG0base"); |
529 | |
530 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
531 | |
532 | TFile::Open(fileNameMC); |
533 | mult->LoadHistograms("Multiplicity"); |
534 | |
535 | const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" }; |
536 | |
537 | TGraph* fitResultsChi2 = new TGraph; |
538 | fitResultsChi2->SetTitle(";Nevents;Chi2"); |
539 | TGraph* fitResultsBayes = new TGraph; |
540 | fitResultsBayes->SetTitle(";Nevents;Chi2"); |
541 | TGraph* fitResultsChi2Limit = new TGraph; |
542 | fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach"); |
543 | TGraph* fitResultsBayesLimit = new TGraph; |
544 | fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach"); |
545 | |
546 | TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600); |
547 | canvas->Divide(5, 2); |
548 | |
549 | Float_t min = 1e10; |
550 | Float_t max = 0; |
551 | |
552 | for (Int_t i=0; i<5; ++i) |
553 | { |
554 | TFile::Open(files[i]); |
555 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
556 | multESD->LoadHistograms("Multiplicity"); |
557 | |
558 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
559 | Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries(); |
560 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(); |
561 | |
562 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3); |
563 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
564 | mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); |
565 | |
566 | Float_t chi2MC = 0; |
567 | Int_t chi2MCLimit = 0; |
568 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
569 | fitResultsChi2->SetPoint(fitResultsChi2->GetN(), nEntries, chi2MC); |
570 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit); |
571 | min = TMath::Min(min, chi2MC); |
572 | max = TMath::Max(max, chi2MC); |
573 | |
574 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result"); |
575 | |
576 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1); |
577 | mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
578 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult"); |
579 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
580 | fitResultsBayes->SetPoint(fitResultsBayes->GetN(), nEntries, chi2MC); |
581 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit); |
582 | |
583 | min = TMath::Min(min, chi2MC); |
584 | max = TMath::Max(max, chi2MC); |
585 | mc->GetXaxis()->SetRangeUser(0, 150); |
586 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
587 | |
588 | // skip errors for now |
589 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
590 | { |
591 | chi2Result->SetBinError(j, 0); |
592 | bayesResult->SetBinError(j, 0); |
593 | } |
594 | |
595 | canvas->cd(i+1); |
596 | mc->SetFillColor(kYellow); |
597 | mc->DrawCopy(); |
598 | chi2Result->SetLineColor(kRed); |
599 | chi2Result->DrawCopy("SAME"); |
600 | bayesResult->SetLineColor(kBlue); |
601 | bayesResult->DrawCopy("SAME"); |
602 | gPad->SetLogy(); |
603 | |
604 | canvas->cd(i+6); |
605 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); |
606 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); |
607 | |
608 | // skip errors for now |
609 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
610 | { |
611 | chi2Result->SetBinError(j, 0); |
612 | bayesResult->SetBinError(j, 0); |
613 | } |
614 | |
615 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); |
616 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
617 | |
618 | chi2Result->DrawCopy(""); |
619 | bayesResult->DrawCopy("SAME"); |
620 | } |
621 | |
622 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); |
623 | canvas->SaveAs(Form("%s.C", canvas->GetName())); |
624 | |
625 | TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400); |
626 | canvas2->Divide(2, 1); |
627 | |
628 | canvas2->cd(1); |
629 | fitResultsChi2->SetMarkerStyle(20); |
630 | fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); |
631 | fitResultsChi2->Draw("AP"); |
632 | |
633 | fitResultsBayes->SetMarkerStyle(3); |
634 | fitResultsBayes->SetMarkerColor(2); |
635 | fitResultsBayes->Draw("P SAME"); |
636 | |
637 | gPad->SetLogy(); |
638 | |
639 | canvas2->cd(2); |
640 | fitResultsChi2Limit->SetMarkerStyle(20); |
641 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); |
642 | fitResultsChi2Limit->Draw("AP"); |
643 | |
644 | fitResultsBayesLimit->SetMarkerStyle(3); |
645 | fitResultsBayesLimit->SetMarkerColor(2); |
646 | fitResultsBayesLimit->Draw("P SAME"); |
647 | |
648 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); |
649 | canvas2->SaveAs(Form("%s.C", canvas2->GetName())); |
650 | } |
651 | |
652 | void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2) |
653 | { |
654 | gSystem->Load("libPWG0base"); |
655 | |
656 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
657 | |
658 | TFile::Open(fileNameMC); |
659 | mult->LoadHistograms("Multiplicity"); |
660 | |
661 | const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" }; |
662 | |
663 | // this one we try to unfold |
664 | TFile::Open(files[0]); |
665 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
666 | multESD->LoadHistograms("Multiplicity"); |
667 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
668 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(); |
669 | |
670 | TGraph* fitResultsChi2 = new TGraph; |
671 | fitResultsChi2->SetTitle(";Input Dist ID;Chi2"); |
672 | TGraph* fitResultsBayes = new TGraph; |
673 | fitResultsBayes->SetTitle(";Input Dist ID;Chi2"); |
674 | TGraph* fitResultsChi2Limit = new TGraph; |
675 | fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach"); |
676 | TGraph* fitResultsBayesLimit = new TGraph; |
677 | fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach"); |
678 | |
679 | TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600); |
680 | canvas->Divide(8, 2); |
681 | |
682 | TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400); |
683 | canvas3->Divide(2, 1); |
684 | |
685 | Float_t min = 1e10; |
686 | Float_t max = 0; |
687 | |
688 | TH1* firstChi = 0; |
689 | TH1* firstBayesian = 0; |
690 | TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond"); |
691 | |
692 | TLegend* legend = new TLegend(0.7, 0.7, 1, 1); |
693 | |
694 | for (Int_t i=0; i<8; ++i) |
695 | { |
696 | if (i == 0) |
697 | { |
698 | startCond = (TH1*) mc->Clone("startCond2"); |
699 | } |
700 | else if (i < 6) |
701 | { |
702 | TFile::Open(files[i-1]); |
703 | AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2"); |
704 | multESD2->LoadHistograms("Multiplicity"); |
705 | startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond"); |
706 | } |
707 | else if (i == 6) |
708 | { |
709 | 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); |
710 | func->SetParNames("scaling", "averagen", "k"); |
711 | func->SetParLimits(0, 1e-3, 1e10); |
712 | func->SetParLimits(1, 0.001, 1000); |
713 | func->SetParLimits(2, 0.001, 1000); |
714 | |
715 | func->SetParameters(1, 10, 2); |
716 | for (Int_t j=2; j<=startCond->GetNbinsX(); j++) |
717 | startCond->SetBinContent(j, func->Eval(j-1)); |
718 | } |
719 | else if (i == 7) |
720 | { |
721 | for (Int_t j=1; j<=startCond->GetNbinsX(); j++) |
722 | startCond->SetBinContent(j, 1); |
723 | } |
724 | |
725 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3); |
726 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond); |
727 | mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); |
728 | |
729 | Float_t chi2MC = 0; |
730 | Int_t chi2MCLimit = 0; |
731 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
732 | fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC); |
733 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit); |
734 | min = TMath::Min(min, chi2MC); |
735 | max = TMath::Max(max, chi2MC); |
736 | |
737 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result"); |
738 | if (!firstChi) |
739 | firstChi = (TH1*) chi2Result->Clone("firstChi"); |
740 | |
741 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1, 30, startCond); |
742 | mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
743 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult"); |
744 | if (!firstBayesian) |
745 | firstBayesian = (TH1*) bayesResult->Clone("firstBayesian"); |
746 | |
747 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
748 | fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC); |
749 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit); |
750 | |
751 | min = TMath::Min(min, chi2MC); |
752 | max = TMath::Max(max, chi2MC); |
753 | mc->GetXaxis()->SetRangeUser(0, 150); |
754 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
755 | |
756 | // skip errors for now |
757 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
758 | { |
759 | chi2Result->SetBinError(j, 0); |
760 | bayesResult->SetBinError(j, 0); |
761 | } |
762 | |
763 | canvas3->cd(1); |
764 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); |
765 | tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); |
766 | tmp->Divide(firstChi); |
767 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
768 | tmp->GetXaxis()->SetRangeUser(0, 200); |
769 | tmp->SetLineColor(i+1); |
770 | legend->AddEntry(tmp, Form("%d", i)); |
771 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
772 | |
773 | canvas3->cd(2); |
774 | tmp = (TH1*) bayesResult->Clone("tmp"); |
775 | tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); |
776 | tmp->Divide(firstBayesian); |
777 | tmp->SetLineColor(i+1); |
778 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
779 | tmp->GetXaxis()->SetRangeUser(0, 200); |
780 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
781 | |
782 | canvas->cd(i+1); |
783 | mc->SetFillColor(kYellow); |
784 | mc->DrawCopy(); |
785 | chi2Result->SetLineColor(kRed); |
786 | chi2Result->DrawCopy("SAME"); |
787 | bayesResult->SetLineColor(kBlue); |
788 | bayesResult->DrawCopy("SAME"); |
789 | gPad->SetLogy(); |
790 | |
791 | canvas->cd(i+9); |
792 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); |
793 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); |
794 | |
795 | // skip errors for now |
796 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
797 | { |
798 | chi2Result->SetBinError(j, 0); |
799 | bayesResult->SetBinError(j, 0); |
800 | } |
801 | |
802 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); |
803 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
804 | |
805 | chi2Result->DrawCopy(""); |
806 | bayesResult->DrawCopy("SAME"); |
807 | } |
808 | |
809 | canvas3->cd(1); |
810 | legend->Draw(); |
811 | |
cfc19dd5 |
812 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); |
dd701109 |
813 | |
814 | TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400); |
815 | canvas2->Divide(2, 1); |
816 | |
817 | canvas2->cd(1); |
818 | fitResultsChi2->SetMarkerStyle(20); |
819 | fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); |
820 | fitResultsChi2->Draw("AP"); |
821 | |
822 | fitResultsBayes->SetMarkerStyle(3); |
823 | fitResultsBayes->SetMarkerColor(2); |
824 | fitResultsBayes->Draw("P SAME"); |
825 | |
826 | gPad->SetLogy(); |
827 | |
828 | canvas2->cd(2); |
829 | fitResultsChi2Limit->SetMarkerStyle(20); |
830 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); |
831 | fitResultsChi2Limit->Draw("AP"); |
832 | |
833 | fitResultsBayesLimit->SetMarkerStyle(3); |
834 | fitResultsBayesLimit->SetMarkerColor(2); |
835 | fitResultsBayesLimit->Draw("P SAME"); |
836 | |
837 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); |
838 | canvas3->SaveAs(Form("%s.gif", canvas3->GetName())); |
cfc19dd5 |
839 | } |
840 | |
841 | void Merge(Int_t n, const char** files, const char* output) |
842 | { |
843 | gSystem->Load("libPWG0base"); |
844 | |
845 | AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n]; |
846 | TList list; |
847 | for (Int_t i=0; i<n; ++i) |
848 | { |
849 | TString name("Multiplicity"); |
850 | if (i > 0) |
851 | name.Form("Multiplicity%d", i); |
852 | |
853 | TFile::Open(files[i]); |
854 | data[i] = new AliMultiplicityCorrection(name, name); |
855 | data[i]->LoadHistograms("Multiplicity"); |
856 | if (i > 0) |
857 | list.Add(data[i]); |
858 | } |
859 | |
860 | data[0]->Merge(&list); |
861 | |
862 | data[0]->DrawHistograms(); |
863 | |
864 | TFile::Open(output, "RECREATE"); |
865 | data[0]->SaveHistograms(); |
866 | gFile->Close(); |
867 | } |
868 | |
869 | void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root") |
870 | { |
871 | gSystem->Load("libPWG0base"); |
872 | |
873 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
874 | |
875 | TFile::Open(fileName); |
876 | mult->LoadHistograms("Multiplicity"); |
877 | |
878 | TF1* func = 0; |
879 | |
880 | if (caseNo >= 4) |
881 | { |
882 | 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); |
883 | func->SetParNames("scaling", "averagen", "k"); |
884 | } |
885 | |
886 | switch (caseNo) |
887 | { |
888 | case 0: func = new TF1("flat", "1"); break; |
889 | case 1: func = new TF1("flat", "501-x"); break; |
890 | case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break; |
891 | case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break; |
dd701109 |
892 | case 4: func->SetParameters(1e7, 10, 2); break; |
893 | case 5: func->SetParameters(1e7, 20, 3); break; |
894 | case 6: func->SetParameters(1e7, 30, 4); break; |
895 | case 7: func->SetParameters(1e7, 70, 2); break; |
896 | case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break; |
cfc19dd5 |
897 | |
898 | default: return; |
899 | } |
900 | |
901 | mult->SetGenMeasFromFunc(func, 2); |
902 | |
dd701109 |
903 | TFile::Open("out.root", "RECREATE"); |
904 | mult->SaveHistograms(); |
905 | |
906 | //mult->ApplyBayesianMethod(2, kFALSE); |
907 | //mult->ApplyMinuitFit(2, kFALSE); |
cfc19dd5 |
908 | //mult->ApplyGaussianMethod(2, kFALSE); |
dd701109 |
909 | mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx); |
910 | } |
911 | |
912 | void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2) |
913 | { |
914 | gSystem->Load("libPWG0base"); |
915 | |
916 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
917 | |
918 | TFile::Open(fileName); |
919 | mult->LoadHistograms("Multiplicity"); |
920 | |
921 | // empty under/overflow bins in x, otherwise Project3D takes them into account |
922 | TH3* corr = mult->GetCorrelation(corrMatrix); |
923 | for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j) |
924 | { |
925 | for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k) |
926 | { |
927 | corr->SetBinContent(0, j, k, 0); |
928 | corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0); |
929 | } |
930 | } |
931 | |
932 | TH2* proj = (TH2*) corr->Project3D("zy"); |
933 | |
934 | // normalize correction for given nPart |
935 | for (Int_t i=1; i<=proj->GetNbinsX(); ++i) |
936 | { |
937 | Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY()); |
938 | if (sum <= 0) |
939 | continue; |
940 | |
941 | for (Int_t j=1; j<=proj->GetNbinsY(); ++j) |
942 | { |
943 | // npart sum to 1 |
944 | proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum); |
945 | proj->SetBinError(i, j, proj->GetBinError(i, j) / sum); |
946 | } |
947 | } |
948 | |
949 | new TCanvas; |
950 | proj->Draw("COLZ"); |
951 | |
952 | TH1* scaling = proj->ProjectionY("scaling", 1, 1); |
953 | scaling->Reset(); |
954 | scaling->SetMarkerStyle(3); |
955 | //scaling->GetXaxis()->SetRangeUser(0, 50); |
956 | TH1* mean = (TH1F*) scaling->Clone("mean"); |
957 | TH1* width = (TH1F*) scaling->Clone("width"); |
958 | |
959 | TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500); |
960 | lognormal->SetParNames("scaling", "mean", "sigma"); |
961 | lognormal->SetParameters(1, 1, 1); |
962 | lognormal->SetParLimits(0, 1, 1); |
963 | lognormal->SetParLimits(1, 0, 100); |
964 | lognormal->SetParLimits(2, 1e-3, 1); |
965 | |
966 | 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); |
967 | nbd->SetParNames("scaling", "averagen", "k"); |
968 | nbd->SetParameters(1, 13, 5); |
969 | nbd->SetParLimits(0, 1, 1); |
970 | nbd->SetParLimits(1, 1, 100); |
971 | nbd->SetParLimits(2, 1, 1e8); |
972 | |
973 | TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50); |
974 | poisson->SetParNames("scaling", "k", "deltax"); |
975 | poisson->SetParameters(1, 1, 0); |
976 | poisson->SetParLimits(0, 0, 10); |
977 | poisson->SetParLimits(1, 0.01, 100); |
978 | poisson->SetParLimits(2, 0, 10); |
979 | |
980 | TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50); |
981 | mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth"); |
982 | mygaus->SetParameters(1, 0, 1, 1, 0, 1); |
983 | mygaus->SetParLimits(2, 1e-5, 10); |
984 | mygaus->SetParLimits(4, 1, 1); |
985 | mygaus->SetParLimits(5, 1e-5, 10); |
986 | |
987 | //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50); |
988 | TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50); |
989 | sqrt->SetParNames("ydelta", "exp", "xdelta"); |
990 | sqrt->SetParameters(0, 0, 1); |
991 | sqrt->SetParLimits(1, 0, 10); |
992 | |
993 | const char* fitWith = "gaus"; |
994 | |
995 | for (Int_t i=1; i<=150; ++i) |
996 | { |
997 | printf("Fitting %d...\n", i); |
998 | |
999 | TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e"); |
1000 | //hist->GetXaxis()->SetRangeUser(0, 50); |
1001 | //lognormal->SetParameter(0, hist->GetMaximum()); |
1002 | hist->Fit(fitWith, "0 M", ""); |
1003 | |
1004 | TF1* func = hist->GetListOfFunctions()->FindObject(fitWith); |
1005 | |
1006 | if (((i-1) % 15 == 0) || ((i % 5 == 0) && i < 30)) |
1007 | { |
1008 | new TCanvas; |
1009 | hist->Draw(); |
1010 | func->Clone()->Draw("SAME"); |
1011 | gPad->SetLogy(); |
1012 | } |
1013 | |
1014 | scaling->Fill(i, func->GetParameter(0)); |
1015 | mean->Fill(i, func->GetParameter(1)); |
1016 | width->Fill(i, func->GetParameter(2)); |
1017 | } |
1018 | |
1019 | TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500); |
1020 | log->SetParameters(0, 1, 1); |
1021 | log->SetParLimits(1, 0, 100); |
1022 | log->SetParLimits(2, 1e-3, 10); |
1023 | |
1024 | TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500); |
1025 | over->SetParameters(0, 1, 0); |
1026 | //over->SetParLimits(0, 0, 100); |
1027 | over->SetParLimits(1, 1e-3, 10); |
1028 | over->SetParLimits(2, 0, 100); |
1029 | |
1030 | c1 = new TCanvas("fitparams", "fitparams", 1200, 400); |
1031 | c1->Divide(3, 1); |
1032 | |
1033 | c1->cd(1); |
1034 | scaling->Draw("P"); |
1035 | |
1036 | //TF1* scalingFit = new TF1("mypol0", "[0]"); |
1037 | TF1* scalingFit = over; |
1038 | scaling->Fit(scalingFit, "", "", 3, 100); |
1039 | |
1040 | c1->cd(2); |
1041 | mean->Draw("P"); |
1042 | |
1043 | //TF1* meanFit = log; |
1044 | TF1* meanFit = new TF1("mypol1", "[0]+[1]*x"); |
1045 | mean->Fit(meanFit, "", "", 3, 100); |
1046 | |
1047 | c1->cd(3); |
1048 | width->Draw("P"); |
1049 | |
1050 | //TF1* widthFit = over; |
1051 | TF1* widthFit = new TF1("mypol2", "[0]+[1]*x+[2]*x*x"); |
1052 | width->Fit(widthFit, "", "", 5, 100); |
1053 | |
1054 | // build new correction matrix |
1055 | TH2* new = (TH2*) proj->Clone("new"); |
1056 | new->Reset(); |
1057 | Float_t x, y; |
1058 | for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) |
1059 | { |
1060 | TF1* func = (TF1*) gROOT->FindObject(fitWith); |
1061 | x = new->GetXaxis()->GetBinCenter(i); |
1062 | //if (i == 1) |
1063 | // x = 0.1; |
1064 | x++; |
1065 | func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); |
1066 | printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); |
1067 | |
1068 | for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) |
1069 | { |
1070 | if (i < 21) |
1071 | { |
1072 | // leave bins 1..20 untouched |
1073 | new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j)); |
1074 | } |
1075 | else |
1076 | { |
1077 | y = new->GetYaxis()->GetBinCenter(j); |
1078 | if (j == 1) |
1079 | y = 0.1; |
1080 | if (func->Eval(y) > 1e-4) |
1081 | new->SetBinContent(i, j, func->Eval(y)); |
1082 | } |
1083 | } |
1084 | } |
1085 | |
1086 | // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0 |
1087 | // we take the values from the old response matrix |
1088 | //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) |
1089 | // new->SetBinContent(i, 1, proj->GetBinContent(i, 1)); |
1090 | |
1091 | //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) |
1092 | // new->SetBinContent(1, j, proj->GetBinContent(1, j)); |
1093 | |
1094 | // normalize correction for given nPart |
1095 | for (Int_t i=1; i<=new->GetNbinsX(); ++i) |
1096 | { |
1097 | Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY()); |
1098 | if (sum <= 0) |
1099 | continue; |
1100 | |
1101 | for (Int_t j=1; j<=new->GetNbinsY(); ++j) |
1102 | { |
1103 | // npart sum to 1 |
1104 | new->SetBinContent(i, j, new->GetBinContent(i, j) / sum); |
1105 | new->SetBinError(i, j, new->GetBinError(i, j) / sum); |
1106 | } |
1107 | } |
1108 | |
1109 | new TCanvas; |
1110 | new->Draw("COLZ"); |
1111 | |
1112 | TH2* diff = (TH2*) new->Clone("diff"); |
1113 | diff->Add(proj, -1); |
1114 | |
1115 | new TCanvas; |
1116 | diff->Draw("COLZ"); |
1117 | diff->SetMinimum(-0.05); |
1118 | diff->SetMaximum(0.05); |
1119 | |
1120 | corr->Reset(); |
1121 | |
1122 | for (Int_t i=1; i<=new->GetNbinsX(); ++i) |
1123 | for (Int_t j=1; j<=new->GetNbinsY(); ++j) |
1124 | corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2, i, j, new->GetBinContent(i, j)); |
1125 | |
1126 | new TCanvas; |
1127 | corr->Project3D("zy")->Draw("COLZ"); |
1128 | |
1129 | TFile::Open("out.root", "RECREATE"); |
1130 | mult->SaveHistograms(); |
1131 | |
1132 | TH1* proj1 = proj->ProjectionY("proj1", 36, 36); |
1133 | TH1* proj2 = new->ProjectionY("proj2", 36, 36); |
1134 | proj2->SetLineColor(2); |
1135 | |
1136 | new TCanvas; |
1137 | proj1->Draw(); |
1138 | proj2->Draw("SAME"); |
cfc19dd5 |
1139 | } |