]>
Commit | Line | Data |
---|---|---|
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 | } |