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) |
447c325d |
13 | connectProof("proof40@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 | |
447c325d |
51 | //Int_t result = chain->Process(selectorName, option); |
03244034 |
52 | Int_t result = executeQuery(chain, &inputList, selectorName, option); |
0ab29cfa |
53 | |
54 | if (result != 0) |
55 | { |
56 | printf("ERROR: Executing process failed with %d.\n", result); |
57 | return; |
58 | } |
0a173978 |
59 | } |
0ab29cfa |
60 | |
0a173978 |
61 | void draw(const char* fileName = "multiplicityMC.root") |
62 | { |
63 | gSystem->Load("libPWG0base"); |
64 | |
65 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
66 | |
67 | TFile::Open(fileName); |
68 | mult->LoadHistograms("Multiplicity"); |
69 | |
70 | mult->DrawHistograms(); |
447c325d |
71 | |
72 | TH2* hist = (TH2*) gROOT->FindObject("fCorrelation3_zy"); |
73 | canvas = new TCanvas("c1", "c1", 600, 500); |
74 | hist->SetStats(kFALSE); |
75 | hist->Draw("COLZ"); |
76 | hist->SetTitle(";true multiplicity in |#eta| < 2;measured multiplicity in |#eta| < 2"); |
77 | hist->GetYaxis()->SetTitleOffset(1.1); |
78 | gPad->SetRightMargin(0.15); |
79 | gPad->SetLogz(); |
80 | |
81 | canvas->SaveAs("Plot_Correlation.pdf"); |
0ab29cfa |
82 | } |
83 | |
cfc19dd5 |
84 | void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0) |
0a173978 |
85 | { |
86 | gSystem->Load("libPWG0base"); |
87 | |
88 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
89 | |
90 | TFile::Open(fileName); |
91 | mult->LoadHistograms("Multiplicity"); |
92 | |
dd701109 |
93 | //mult->ApplyLaszloMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx); |
94 | |
95 | //return; |
96 | |
97 | |
447c325d |
98 | /*mult->ApplyBayesianMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx); |
dd701109 |
99 | mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY()); |
447c325d |
100 | return;*/ |
dd701109 |
101 | |
102 | TStopwatch timer; |
cfc19dd5 |
103 | |
dd701109 |
104 | timer.Start(); |
105 | |
447c325d |
106 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.2); |
107 | mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); |
dd701109 |
108 | mult->DrawComparison("MinuitChi2", hist, kFALSE, kFALSE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY()); |
109 | |
110 | timer.Stop(); |
111 | timer.Print(); |
112 | |
113 | return 0; |
cfc19dd5 |
114 | |
115 | //mult->ApplyGaussianMethod(hist, kFALSE); |
9ca6be09 |
116 | |
dd701109 |
117 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); |
118 | mult->ApplyNBDFit(hist, kFALSE); |
119 | mult->DrawComparison("NBDChi2Fit", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY()); |
120 | |
9ca6be09 |
121 | return mult; |
122 | } |
123 | |
447c325d |
124 | void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Bool_t chi2 = kTRUE, Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE) |
9ca6be09 |
125 | { |
126 | gSystem->Load("libPWG0base"); |
127 | |
128 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
129 | |
130 | TFile::Open(fileNameMC); |
131 | mult->LoadHistograms("Multiplicity"); |
132 | |
133 | TFile::Open(fileNameESD); |
cfc19dd5 |
134 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); |
447c325d |
135 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID))); |
136 | //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID)); |
9ca6be09 |
137 | |
cfc19dd5 |
138 | mult->SetMultiplicityESD(histID, hist); |
9ca6be09 |
139 | |
447c325d |
140 | // small hack to get around charge conservation for full phase space ;-) |
141 | if (fullPhaseSpace) |
142 | { |
143 | TH1* corr = mult->GetCorrelation(histID + 4); |
dd701109 |
144 | |
447c325d |
145 | for (Int_t i=2; i<=corr->GetNbinsX(); i+=2) |
146 | for (Int_t j=1; j<=corr->GetNbinsY(); ++j) |
147 | { |
148 | corr->SetBinContent(i, j, corr->GetBinContent(i-1, j)); |
149 | corr->SetBinError(i, j, corr->GetBinError(i-1, j)); |
150 | } |
151 | } |
dd701109 |
152 | |
447c325d |
153 | /*mult->SetMultiplicityVtx(histID, hist2); |
154 | mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
155 | return;*/ |
cfc19dd5 |
156 | |
447c325d |
157 | if (chi2) |
158 | { |
159 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4); |
160 | mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); //, hist2->ProjectionY("mymchist")); |
161 | mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist")); |
162 | } |
163 | else |
dd701109 |
164 | { |
447c325d |
165 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100); |
166 | mult->DrawComparison("Bayesian", histID, kFALSE, kTRUE, hist2->ProjectionY("mymchist2")); |
dd701109 |
167 | } |
168 | |
169 | //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e7); |
170 | //mult->ApplyMinuitFit(histID, kFALSE); |
171 | //mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY()); |
172 | |
447c325d |
173 | } |
174 | |
175 | void* fit2Step(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE) |
176 | { |
177 | gSystem->Load("libPWG0base"); |
178 | |
179 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
180 | |
181 | TFile::Open(fileNameMC); |
182 | mult->LoadHistograms("Multiplicity"); |
183 | |
184 | TFile::Open(fileNameESD); |
185 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); |
186 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID))); |
187 | //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID)); |
188 | |
189 | mult->SetMultiplicityESD(histID, hist); |
190 | |
191 | // small hack to get around charge conservation for full phase space ;-) |
192 | if (fullPhaseSpace) |
193 | { |
194 | TH1* corr = mult->GetCorrelation(histID + 4); |
195 | |
196 | for (Int_t i=2; i<=corr->GetNbinsX(); i+=2) |
197 | for (Int_t j=1; j<=corr->GetNbinsY(); ++j) |
198 | { |
199 | corr->SetBinContent(i, j, corr->GetBinContent(i-1, j)); |
200 | corr->SetBinError(i, j, corr->GetBinError(i-1, j)); |
201 | } |
202 | } |
203 | |
204 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
205 | mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); |
206 | mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist")); |
207 | |
208 | TH1* result = (TH1*) mult->GetMultiplicityESDCorrected((fullPhaseSpace) ? 4 : histID))->Clone("firstresult"); |
209 | |
210 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 100000); |
211 | mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE, result); |
212 | mult->DrawComparison("MinuitChi2_Step2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist")); |
0a173978 |
213 | |
214 | return mult; |
215 | } |
cfc19dd5 |
216 | |
217 | const char* GetRegName(Int_t type) |
218 | { |
219 | switch (type) |
220 | { |
221 | case AliMultiplicityCorrection::kNone: return "None"; break; |
222 | case AliMultiplicityCorrection::kPol0: return "Pol0"; break; |
223 | case AliMultiplicityCorrection::kPol1: return "Pol1"; break; |
224 | case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break; |
225 | case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break; |
dd701109 |
226 | case AliMultiplicityCorrection::kTest : return "Test"; break; |
227 | } |
228 | return 0; |
229 | } |
230 | |
231 | const char* GetEventTypeName(Int_t type) |
232 | { |
233 | switch (type) |
234 | { |
235 | case AliMultiplicityCorrection::kTrVtx: return "trigger, vertex"; break; |
236 | case AliMultiplicityCorrection::kMB: return "minimum bias"; break; |
237 | case AliMultiplicityCorrection::kINEL: return "inelastic"; break; |
cfc19dd5 |
238 | } |
239 | return 0; |
240 | } |
241 | |
447c325d |
242 | void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3) |
cfc19dd5 |
243 | { |
dd701109 |
244 | gSystem->mkdir(targetDir); |
245 | |
246 | gSystem->Load("libPWG0base"); |
247 | |
248 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
249 | TFile::Open(fileNameMC); |
250 | mult->LoadHistograms("Multiplicity"); |
251 | |
252 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
253 | TFile::Open(fileNameESD); |
254 | multESD->LoadHistograms("Multiplicity"); |
255 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
256 | |
257 | TCanvas* canvas = new TCanvas("EvaluateBayesianMethod", "EvaluateBayesianMethod", 800, 600); |
258 | TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); |
259 | legend->SetFillColor(0); |
260 | |
261 | Float_t min = 1e10; |
262 | Float_t max = 0; |
263 | |
264 | TGraph* first = 0; |
265 | Int_t count = 0; // just to order the saved images... |
266 | |
267 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) |
268 | { |
269 | TGraph* fitResultsMC = new TGraph; |
270 | fitResultsMC->SetTitle(";Weight Parameter"); |
271 | TGraph* fitResultsRes = new TGraph; |
272 | fitResultsRes->SetTitle(";Weight Parameter"); |
273 | |
274 | fitResultsMC->SetFillColor(0); |
275 | fitResultsRes->SetFillColor(0); |
276 | fitResultsMC->SetMarkerStyle(20+type); |
277 | fitResultsRes->SetMarkerStyle(24+type); |
278 | fitResultsRes->SetMarkerColor(kRed); |
279 | fitResultsRes->SetLineColor(kRed); |
280 | |
281 | legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type))); |
282 | legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type))); |
283 | |
447c325d |
284 | for (Float_t weight = 0; weight < 1.01; weight += 0.1) |
dd701109 |
285 | { |
286 | Float_t chi2MC = 0; |
287 | Float_t residuals = 0; |
288 | |
289 | mult->ApplyBayesianMethod(histID, kFALSE, type, weight); |
290 | mult->DrawComparison(Form("%s/Bayesian_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); |
291 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
292 | |
293 | fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC); |
294 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); |
295 | |
296 | min = TMath::Min(min, TMath::Min(chi2MC, residuals)); |
297 | max = TMath::Max(max, TMath::Max(chi2MC, residuals)); |
298 | } |
299 | |
300 | fitResultsMC->Print(); |
301 | fitResultsRes->Print(); |
302 | |
303 | canvas->cd(); |
304 | fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME")); |
305 | fitResultsRes->Draw("SAME CP"); |
306 | |
307 | if (first == 0) |
308 | first = fitResultsMC; |
309 | } |
310 | |
311 | gPad->SetLogy(); |
312 | printf("min = %f, max = %f\n", min, max); |
313 | if (min <= 0) |
314 | min = 1e-5; |
315 | first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5); |
316 | |
317 | legend->Draw(); |
318 | |
319 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
320 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); |
321 | } |
322 | |
323 | void EvaluateBayesianMethodIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) |
324 | { |
325 | gSystem->mkdir(targetDir); |
326 | |
327 | gSystem->Load("libPWG0base"); |
328 | |
329 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
330 | TFile::Open(fileNameMC); |
331 | mult->LoadHistograms("Multiplicity"); |
332 | |
333 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
334 | TFile::Open(fileNameESD); |
335 | multESD->LoadHistograms("Multiplicity"); |
336 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
337 | |
338 | TCanvas* canvas = new TCanvas("EvaluateBayesianMethodIterations", "EvaluateBayesianMethodIterations", 800, 600); |
339 | TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); |
340 | legend->SetFillColor(0); |
341 | |
342 | Float_t min = 1e10; |
343 | Float_t max = 0; |
344 | |
345 | TGraph* first = 0; |
346 | Int_t count = 0; // just to order the saved images... |
347 | |
348 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) |
349 | { |
350 | TGraph* fitResultsMC = new TGraph; |
351 | fitResultsMC->SetTitle(";Iterations"); |
352 | TGraph* fitResultsRes = new TGraph; |
353 | fitResultsRes->SetTitle(";Iterations"); |
354 | |
355 | fitResultsMC->SetFillColor(0); |
356 | fitResultsRes->SetFillColor(0); |
357 | fitResultsMC->SetMarkerStyle(20+type); |
358 | fitResultsRes->SetMarkerStyle(24+type); |
359 | fitResultsRes->SetMarkerColor(kRed); |
360 | fitResultsRes->SetLineColor(kRed); |
361 | |
362 | legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type))); |
363 | legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type))); |
364 | |
365 | for (Int_t iter = 5; iter <= 50; iter += 5) |
366 | { |
367 | Float_t chi2MC = 0; |
368 | Float_t residuals = 0; |
369 | |
370 | mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1, iter); |
371 | mult->DrawComparison(Form("%s/BayesianIter_%02d_%d_%d", targetDir, count++, type, iter), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); |
372 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
373 | |
374 | fitResultsMC->SetPoint(fitResultsMC->GetN(), iter, chi2MC); |
375 | fitResultsRes->SetPoint(fitResultsRes->GetN(), iter, residuals); |
376 | |
377 | min = TMath::Min(min, TMath::Min(chi2MC, residuals)); |
378 | max = TMath::Max(max, TMath::Max(chi2MC, residuals)); |
379 | } |
380 | |
381 | fitResultsMC->Print(); |
382 | fitResultsRes->Print(); |
383 | |
384 | canvas->cd(); |
385 | fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME")); |
386 | fitResultsRes->Draw("SAME CP"); |
387 | |
388 | if (first == 0) |
389 | first = fitResultsMC; |
390 | } |
391 | |
392 | gPad->SetLogy(); |
393 | printf("min = %f, max = %f\n", min, max); |
394 | if (min <= 0) |
395 | min = 1e-5; |
396 | first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5); |
397 | |
398 | legend->Draw(); |
399 | |
400 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
401 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); |
402 | } |
403 | |
447c325d |
404 | void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3) |
405 | { |
406 | gSystem->mkdir(targetDir); |
407 | |
408 | gSystem->Load("libPWG0base"); |
409 | |
410 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
411 | TFile::Open(fileNameMC); |
412 | mult->LoadHistograms("Multiplicity"); |
413 | |
414 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
415 | TFile::Open(fileNameESD); |
416 | multESD->LoadHistograms("Multiplicity"); |
417 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
418 | |
419 | Int_t count = 0; // just to order the saved images... |
420 | |
421 | TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE"); |
422 | |
423 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) |
424 | { |
425 | TString tmp; |
426 | tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type)); |
427 | |
428 | TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600); |
429 | |
430 | for (Int_t i = 1; i <= 7; i++) |
431 | { |
432 | Int_t iter = i * 20; |
433 | TGraph* fitResultsMC = new TGraph; |
434 | fitResultsMC->SetTitle(Form("%d iterations", iter)); |
435 | fitResultsMC->GetXaxis()->SetTitle("smoothing parameter"); |
436 | fitResultsMC->GetYaxis()->SetTitle("P_{1} (2 <= t <= 150)"); |
437 | fitResultsMC->SetName(Form("%s_MC_%d", tmp.Data(), i)); |
438 | TGraph* fitResultsRes = new TGraph; |
439 | fitResultsRes->SetTitle(Form("%d iterations", iter)); |
440 | fitResultsRes->SetName(Form("%s_Res_%d", tmp.Data(), i)); |
441 | fitResultsRes->GetXaxis()->SetTitle("smoothing parameter"); |
442 | fitResultsRes->GetYaxis()->SetTitle("P_{2}"); |
443 | |
444 | fitResultsMC->SetFillColor(0); |
445 | fitResultsRes->SetFillColor(0); |
446 | fitResultsMC->SetMarkerStyle(19+i); |
447 | fitResultsRes->SetMarkerStyle(19+i); |
448 | fitResultsRes->SetMarkerColor(kRed); |
449 | fitResultsRes->SetLineColor(kRed); |
450 | |
451 | for (Float_t weight = 0.0; weight < 1.01; weight += 0.2) |
452 | { |
453 | Float_t chi2MC = 0; |
454 | Float_t residuals = 0; |
455 | |
456 | mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter); |
457 | mult->DrawComparison(Form("%s/BayesianIterSmooth_%03d_%d_%d_%f", targetDir, count++, type, iter, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); |
458 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
459 | |
460 | fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC); |
461 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); |
462 | } |
463 | |
464 | fitResultsMC->Write(); |
465 | fitResultsRes->Write(); |
466 | } |
467 | } |
468 | } |
469 | |
470 | void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE) |
471 | { |
472 | gSystem->Load("libPWG0base"); |
473 | |
474 | TString name; |
475 | TFile* graphFile = 0; |
476 | if (type == -1) |
477 | { |
478 | name = "EvaluateChi2Method"; |
479 | graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir)); |
480 | } |
481 | else |
482 | { |
483 | name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type)); |
484 | graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir)); |
485 | } |
486 | |
487 | TCanvas* canvas = new TCanvas(name, name, 800, 500); |
488 | if (type == -1) |
489 | canvas->SetLogx(); |
490 | canvas->SetLogy(); |
491 | |
492 | TLegend* legend = new TLegend(0.6, 0.75, 0.88, 0.88); |
493 | legend->SetFillColor(0); |
494 | |
495 | Int_t count = 1; |
496 | |
497 | Float_t xMin = 1e20; |
498 | Float_t xMax = 0; |
499 | |
500 | Float_t yMin = 1e20; |
501 | Float_t yMax = 0; |
502 | |
503 | TString xaxis, yaxis; |
504 | |
505 | while (1) |
506 | { |
507 | TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count)); |
508 | TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count)); |
509 | |
510 | if (!mc || !res) |
511 | break; |
512 | |
513 | xaxis = mc->GetXaxis()->GetTitle(); |
514 | yaxis = mc->GetYaxis()->GetTitle(); |
515 | |
516 | mc->Print(); |
517 | res->Print(); |
518 | |
519 | count++; |
520 | |
521 | if (plotRes) |
522 | { |
523 | xMin = TMath::Min(TMath::Min(xMin, mc->GetXaxis()->GetXmin()), res->GetXaxis()->GetXmin()); |
524 | yMin = TMath::Min(TMath::Min(yMin, mc->GetYaxis()->GetXmin()), res->GetYaxis()->GetXmin()); |
525 | |
526 | xMax = TMath::Max(TMath::Max(xMax, mc->GetXaxis()->GetXmax()), res->GetXaxis()->GetXmax()); |
527 | yMax = TMath::Max(TMath::Max(yMax, mc->GetYaxis()->GetXmax()), res->GetYaxis()->GetXmax()); |
528 | } |
529 | else |
530 | { |
531 | xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin()); |
532 | yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin()); |
533 | |
534 | xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax()); |
535 | yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax()); |
536 | } |
537 | } |
538 | |
539 | if (type >= 0) |
540 | { |
541 | xaxis = "smoothing parameter"; |
542 | } |
543 | else if (type == -1) |
544 | { |
545 | xaxis = "weight parameter"; |
546 | //yaxis = "P_{3} (2 <= t <= 150)"; |
547 | } |
548 | yaxis = "P_{1} (2 <= t <= 150)"; |
549 | |
550 | printf("%f %f %f %f\n", xMin, xMax, yMin, yMax); |
551 | |
552 | TGraph* dummy = new TGraph; |
553 | dummy->SetPoint(0, xMin, yMin); |
554 | dummy->SetPoint(1, xMax, yMax); |
555 | dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data())); |
556 | |
557 | dummy->SetMarkerColor(0); |
558 | dummy->Draw("AP"); |
559 | |
560 | count = 1; |
561 | |
562 | while (1) |
563 | { |
564 | TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count)); |
565 | TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count)); |
566 | |
567 | printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res); |
568 | |
569 | if (!mc || !res) |
570 | break; |
571 | |
572 | printf("Loaded %d sucessful.\n", count); |
573 | |
574 | if (type == -1) |
575 | { |
576 | legend->AddEntry(mc, Form("Eq. (%d)", count+9)); |
577 | } |
578 | else |
579 | legend->AddEntry(mc); |
580 | |
581 | mc->Draw("SAME PC"); |
582 | |
583 | if (plotRes) |
584 | { |
585 | legend->AddEntry(res); |
586 | res->Draw("SAME PC"); |
587 | } |
588 | |
589 | count++; |
590 | } |
591 | |
592 | legend->Draw(); |
593 | |
594 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
595 | canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName())); |
596 | } |
597 | |
598 | void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3) |
dd701109 |
599 | { |
600 | gSystem->mkdir(targetDir); |
601 | |
cfc19dd5 |
602 | gSystem->Load("libPWG0base"); |
603 | |
604 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
605 | |
606 | TFile::Open(fileNameMC); |
607 | mult->LoadHistograms("Multiplicity"); |
608 | |
609 | TFile::Open(fileNameESD); |
610 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); |
611 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID)); |
612 | |
613 | mult->SetMultiplicityESD(histID, hist); |
614 | |
dd701109 |
615 | Int_t count = 0; // just to order the saved images... |
616 | |
447c325d |
617 | TGraph* fitResultsMC = 0; |
618 | TGraph* fitResultsRes = 0; |
619 | |
620 | TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE"); |
cfc19dd5 |
621 | |
622 | for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type) |
447c325d |
623 | // for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type) |
cfc19dd5 |
624 | { |
447c325d |
625 | fitResultsMC = new TGraph; |
626 | fitResultsMC->SetTitle(Form("%s", GetRegName(type))); |
627 | fitResultsMC->GetXaxis()->SetTitle("Weight Parameter"); |
628 | fitResultsMC->SetName(Form("EvaluateChi2Method_MC_%d", type)); |
629 | fitResultsRes = new TGraph; |
630 | fitResultsRes->SetTitle(Form("%s residual chi2", GetRegName(type))); |
631 | fitResultsRes->SetName(Form("EvaluateChi2Method_Res_%d", type)); |
632 | fitResultsRes->GetXaxis()->SetTitle("Weight Parameter"); |
cfc19dd5 |
633 | |
634 | fitResultsMC->SetFillColor(0); |
635 | fitResultsRes->SetFillColor(0); |
636 | fitResultsMC->SetMarkerStyle(19+type); |
dd701109 |
637 | fitResultsRes->SetMarkerStyle(23+type); |
cfc19dd5 |
638 | fitResultsRes->SetMarkerColor(kRed); |
639 | fitResultsRes->SetLineColor(kRed); |
640 | |
447c325d |
641 | for (Int_t i=0; i<5; ++i) |
642 | { |
643 | Float_t weight = TMath::Power(TMath::Sqrt(10), i+6); |
644 | //Float_t weight = TMath::Power(10, i+2); |
cfc19dd5 |
645 | |
447c325d |
646 | //if (type == AliMultiplicityCorrection::kEntropy) weight = 1e4 * (i+1) * 1.5; |
cfc19dd5 |
647 | |
cfc19dd5 |
648 | Float_t chi2MC = 0; |
649 | Float_t residuals = 0; |
447c325d |
650 | Float_t chi2Limit = 0; |
651 | |
652 | TString runName; |
653 | runName.Form("MinuitChi2_%02d_%d_%f", count++, type, weight); |
cfc19dd5 |
654 | |
655 | mult->SetRegularizationParameters(type, weight); |
447c325d |
656 | Int_t status = mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
657 | mult->DrawComparison(Form("%s/%s", targetDir, runName.Data()), histID, kFALSE, kTRUE, hist2->ProjectionY()); |
658 | if (status != 0) |
659 | { |
660 | printf("MINUIT did not succeed. Skipping...\n"); |
661 | continue; |
662 | } |
663 | |
dd701109 |
664 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
447c325d |
665 | TH1* result = mult->GetMultiplicityESDCorrected(histID); |
666 | result->SetName(runName); |
667 | result->Write(); |
cfc19dd5 |
668 | |
669 | fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC); |
670 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); |
cfc19dd5 |
671 | } |
672 | |
447c325d |
673 | graphFile->cd(); |
674 | fitResultsMC->Write(); |
675 | fitResultsRes->Write(); |
cfc19dd5 |
676 | } |
677 | |
447c325d |
678 | graphFile->Close(); |
dd701109 |
679 | } |
680 | |
681 | void EvaluateChi2MethodAll() |
682 | { |
683 | EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); |
684 | EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); |
685 | EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); |
686 | EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); |
687 | EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); |
688 | } |
689 | |
690 | void EvaluateBayesianMethodAll() |
691 | { |
692 | EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); |
693 | EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); |
694 | EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); |
695 | EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); |
696 | EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); |
697 | } |
698 | |
447c325d |
699 | void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3) |
dd701109 |
700 | { |
701 | gSystem->mkdir(targetDir); |
702 | |
703 | gSystem->Load("libPWG0base"); |
704 | |
705 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
706 | |
707 | TFile::Open(fileNameMC); |
708 | mult->LoadHistograms("Multiplicity"); |
709 | |
710 | TFile::Open(fileNameESD); |
711 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
712 | multESD->LoadHistograms("Multiplicity"); |
713 | |
714 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
715 | |
447c325d |
716 | TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 1200); |
717 | canvas->Divide(3, 3); |
dd701109 |
718 | |
719 | Int_t count = 0; |
720 | |
447c325d |
721 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type) |
dd701109 |
722 | { |
447c325d |
723 | TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY("mymc"); |
724 | mc->Sumw2(); |
dd701109 |
725 | |
447c325d |
726 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
dd701109 |
727 | mult->ApplyMinuitFit(histID, kFALSE, type); |
728 | mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc); |
729 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result"); |
730 | |
731 | mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1); |
732 | mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc); |
733 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult"); |
734 | |
735 | mc->GetXaxis()->SetRangeUser(0, 150); |
736 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
737 | |
447c325d |
738 | /* // skip errors for now |
dd701109 |
739 | for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i) |
740 | { |
741 | chi2Result->SetBinError(i, 0); |
742 | bayesResult->SetBinError(i, 0); |
447c325d |
743 | }*/ |
dd701109 |
744 | |
745 | canvas->cd(++count); |
746 | mc->SetFillColor(kYellow); |
747 | mc->DrawCopy(); |
748 | chi2Result->SetLineColor(kRed); |
749 | chi2Result->DrawCopy("SAME"); |
750 | bayesResult->SetLineColor(kBlue); |
751 | bayesResult->DrawCopy("SAME"); |
752 | gPad->SetLogy(); |
753 | |
754 | canvas->cd(count + 3); |
447c325d |
755 | chi2ResultRatio = (TH1*) chi2Result->Clone("chi2ResultRatio"); |
756 | bayesResultRatio = (TH1*) bayesResult->Clone("bayesResultRatio"); |
757 | chi2ResultRatio->Divide(chi2Result, mc, 1, 1, ""); |
758 | bayesResultRatio->Divide(bayesResult, mc, 1, 1, ""); |
dd701109 |
759 | |
447c325d |
760 | chi2ResultRatio->GetYaxis()->SetRangeUser(0.5, 1.5); |
dd701109 |
761 | |
447c325d |
762 | chi2ResultRatio->DrawCopy("HIST"); |
763 | bayesResultRatio->DrawCopy("SAME HIST"); |
dd701109 |
764 | |
447c325d |
765 | canvas->cd(count + 6); |
766 | chi2Result->Divide(chi2Result, bayesResult, 1, 1, ""); |
767 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
768 | chi2Result->DrawCopy("HIST"); |
dd701109 |
769 | } |
770 | |
771 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); |
772 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); |
773 | } |
774 | |
447c325d |
775 | void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3) |
dd701109 |
776 | { |
777 | gSystem->Load("libPWG0base"); |
778 | |
779 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
780 | |
781 | TFile::Open(fileNameMC); |
782 | mult->LoadHistograms("Multiplicity"); |
783 | |
447c325d |
784 | const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" }; |
dd701109 |
785 | |
447c325d |
786 | TGraph* fitResultsChi2 = new TGraph; fitResultsChi2->SetTitle(";Nevents;Chi2"); |
787 | TGraph* fitResultsBayes = new TGraph; fitResultsBayes->SetTitle(";Nevents;Chi2"); |
788 | TGraph* fitResultsChi2Limit = new TGraph; fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach"); |
789 | TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach"); |
790 | |
791 | fitResultsChi2->SetName("fitResultsChi2"); |
792 | fitResultsBayes->SetName("fitResultsBayes"); |
793 | fitResultsChi2Limit->SetName("fitResultsChi2Limit"); |
794 | fitResultsBayesLimit->SetName("fitResultsBayesLimit"); |
dd701109 |
795 | |
796 | TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600); |
797 | canvas->Divide(5, 2); |
798 | |
799 | Float_t min = 1e10; |
800 | Float_t max = 0; |
801 | |
447c325d |
802 | TFile* file = TFile::Open("StatisticsPlot.root", "RECREATE"); |
803 | file->Close(); |
804 | |
dd701109 |
805 | for (Int_t i=0; i<5; ++i) |
806 | { |
807 | TFile::Open(files[i]); |
808 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
809 | multESD->LoadHistograms("Multiplicity"); |
810 | |
811 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
812 | Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries(); |
447c325d |
813 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i)); |
dd701109 |
814 | |
447c325d |
815 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
dd701109 |
816 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
817 | mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); |
818 | |
819 | Float_t chi2MC = 0; |
820 | Int_t chi2MCLimit = 0; |
821 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
822 | fitResultsChi2->SetPoint(fitResultsChi2->GetN(), nEntries, chi2MC); |
823 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit); |
824 | min = TMath::Min(min, chi2MC); |
825 | max = TMath::Max(max, chi2MC); |
826 | |
447c325d |
827 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i)); |
dd701109 |
828 | |
447c325d |
829 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.07, 10); |
dd701109 |
830 | mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
447c325d |
831 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i)); |
dd701109 |
832 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
833 | fitResultsBayes->SetPoint(fitResultsBayes->GetN(), nEntries, chi2MC); |
834 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit); |
835 | |
836 | min = TMath::Min(min, chi2MC); |
837 | max = TMath::Max(max, chi2MC); |
838 | mc->GetXaxis()->SetRangeUser(0, 150); |
839 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
840 | |
841 | // skip errors for now |
842 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
843 | { |
844 | chi2Result->SetBinError(j, 0); |
845 | bayesResult->SetBinError(j, 0); |
846 | } |
847 | |
848 | canvas->cd(i+1); |
849 | mc->SetFillColor(kYellow); |
850 | mc->DrawCopy(); |
851 | chi2Result->SetLineColor(kRed); |
852 | chi2Result->DrawCopy("SAME"); |
853 | bayesResult->SetLineColor(kBlue); |
854 | bayesResult->DrawCopy("SAME"); |
855 | gPad->SetLogy(); |
856 | |
857 | canvas->cd(i+6); |
858 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); |
859 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); |
860 | |
861 | // skip errors for now |
862 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
863 | { |
864 | chi2Result->SetBinError(j, 0); |
865 | bayesResult->SetBinError(j, 0); |
866 | } |
867 | |
868 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); |
869 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
870 | |
871 | chi2Result->DrawCopy(""); |
872 | bayesResult->DrawCopy("SAME"); |
447c325d |
873 | |
874 | TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE"); |
875 | mc->Write(); |
876 | chi2Result->Write(); |
877 | bayesResult->Write(); |
878 | file->Close(); |
dd701109 |
879 | } |
880 | |
881 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); |
882 | canvas->SaveAs(Form("%s.C", canvas->GetName())); |
883 | |
884 | TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400); |
885 | canvas2->Divide(2, 1); |
886 | |
887 | canvas2->cd(1); |
888 | fitResultsChi2->SetMarkerStyle(20); |
889 | fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); |
890 | fitResultsChi2->Draw("AP"); |
891 | |
892 | fitResultsBayes->SetMarkerStyle(3); |
893 | fitResultsBayes->SetMarkerColor(2); |
894 | fitResultsBayes->Draw("P SAME"); |
895 | |
896 | gPad->SetLogy(); |
897 | |
898 | canvas2->cd(2); |
899 | fitResultsChi2Limit->SetMarkerStyle(20); |
900 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); |
901 | fitResultsChi2Limit->Draw("AP"); |
902 | |
903 | fitResultsBayesLimit->SetMarkerStyle(3); |
904 | fitResultsBayesLimit->SetMarkerColor(2); |
905 | fitResultsBayesLimit->Draw("P SAME"); |
906 | |
907 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); |
908 | canvas2->SaveAs(Form("%s.C", canvas2->GetName())); |
447c325d |
909 | |
910 | TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE"); |
911 | fitResultsChi2->Write(); |
912 | fitResultsBayes->Write(); |
913 | fitResultsChi2Limit->Write(); |
914 | fitResultsBayesLimit->Write(); |
915 | file->Close(); |
dd701109 |
916 | } |
917 | |
447c325d |
918 | void StartingConditions(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3) |
dd701109 |
919 | { |
920 | gSystem->Load("libPWG0base"); |
921 | |
922 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
923 | |
924 | TFile::Open(fileNameMC); |
925 | mult->LoadHistograms("Multiplicity"); |
926 | |
447c325d |
927 | const char* files[] = { "multiplicityMC_1M_3.root", "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root" } |
dd701109 |
928 | |
929 | // this one we try to unfold |
930 | TFile::Open(files[0]); |
931 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); |
932 | multESD->LoadHistograms("Multiplicity"); |
933 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
447c325d |
934 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY("mc"); |
dd701109 |
935 | |
936 | TGraph* fitResultsChi2 = new TGraph; |
937 | fitResultsChi2->SetTitle(";Input Dist ID;Chi2"); |
938 | TGraph* fitResultsBayes = new TGraph; |
939 | fitResultsBayes->SetTitle(";Input Dist ID;Chi2"); |
940 | TGraph* fitResultsChi2Limit = new TGraph; |
941 | fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach"); |
942 | TGraph* fitResultsBayesLimit = new TGraph; |
943 | fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach"); |
944 | |
945 | TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600); |
946 | canvas->Divide(8, 2); |
947 | |
948 | TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400); |
949 | canvas3->Divide(2, 1); |
950 | |
951 | Float_t min = 1e10; |
952 | Float_t max = 0; |
953 | |
954 | TH1* firstChi = 0; |
955 | TH1* firstBayesian = 0; |
956 | TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond"); |
957 | |
958 | TLegend* legend = new TLegend(0.7, 0.7, 1, 1); |
959 | |
447c325d |
960 | TFile* file = TFile::Open("StartingConditions.root", "RECREATE"); |
961 | mc->Write(); |
962 | file->Close(); |
963 | |
dd701109 |
964 | for (Int_t i=0; i<8; ++i) |
965 | { |
966 | if (i == 0) |
967 | { |
968 | startCond = (TH1*) mc->Clone("startCond2"); |
969 | } |
970 | else if (i < 6) |
971 | { |
972 | TFile::Open(files[i-1]); |
973 | AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2"); |
974 | multESD2->LoadHistograms("Multiplicity"); |
975 | startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond"); |
976 | } |
977 | else if (i == 6) |
978 | { |
979 | 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); |
980 | func->SetParNames("scaling", "averagen", "k"); |
981 | func->SetParLimits(0, 1e-3, 1e10); |
982 | func->SetParLimits(1, 0.001, 1000); |
983 | func->SetParLimits(2, 0.001, 1000); |
984 | |
985 | func->SetParameters(1, 10, 2); |
986 | for (Int_t j=2; j<=startCond->GetNbinsX(); j++) |
987 | startCond->SetBinContent(j, func->Eval(j-1)); |
988 | } |
989 | else if (i == 7) |
990 | { |
991 | for (Int_t j=1; j<=startCond->GetNbinsX(); j++) |
992 | startCond->SetBinContent(j, 1); |
993 | } |
994 | |
447c325d |
995 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
dd701109 |
996 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond); |
997 | mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); |
998 | |
999 | Float_t chi2MC = 0; |
1000 | Int_t chi2MCLimit = 0; |
1001 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
1002 | fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC); |
1003 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit); |
1004 | min = TMath::Min(min, chi2MC); |
1005 | max = TMath::Max(max, chi2MC); |
1006 | |
447c325d |
1007 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i)); |
dd701109 |
1008 | if (!firstChi) |
1009 | firstChi = (TH1*) chi2Result->Clone("firstChi"); |
1010 | |
447c325d |
1011 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, startCond); |
dd701109 |
1012 | mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
447c325d |
1013 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i)); |
dd701109 |
1014 | if (!firstBayesian) |
1015 | firstBayesian = (TH1*) bayesResult->Clone("firstBayesian"); |
1016 | |
1017 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
1018 | fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC); |
1019 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit); |
1020 | |
447c325d |
1021 | TFile* file = TFile::Open("StartingConditions.root", "UPDATE"); |
1022 | chi2Result->Write(); |
1023 | bayesResult->Write(); |
1024 | file->Close(); |
1025 | |
dd701109 |
1026 | min = TMath::Min(min, chi2MC); |
1027 | max = TMath::Max(max, chi2MC); |
1028 | mc->GetXaxis()->SetRangeUser(0, 150); |
1029 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
1030 | |
1031 | // skip errors for now |
1032 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
1033 | { |
1034 | chi2Result->SetBinError(j, 0); |
1035 | bayesResult->SetBinError(j, 0); |
1036 | } |
1037 | |
1038 | canvas3->cd(1); |
1039 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); |
1040 | tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); |
1041 | tmp->Divide(firstChi); |
1042 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1043 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1044 | tmp->SetLineColor(i+1); |
1045 | legend->AddEntry(tmp, Form("%d", i)); |
1046 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1047 | |
1048 | canvas3->cd(2); |
1049 | tmp = (TH1*) bayesResult->Clone("tmp"); |
1050 | tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); |
1051 | tmp->Divide(firstBayesian); |
1052 | tmp->SetLineColor(i+1); |
1053 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1054 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1055 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1056 | |
1057 | canvas->cd(i+1); |
1058 | mc->SetFillColor(kYellow); |
1059 | mc->DrawCopy(); |
1060 | chi2Result->SetLineColor(kRed); |
1061 | chi2Result->DrawCopy("SAME"); |
1062 | bayesResult->SetLineColor(kBlue); |
1063 | bayesResult->DrawCopy("SAME"); |
1064 | gPad->SetLogy(); |
1065 | |
1066 | canvas->cd(i+9); |
1067 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); |
1068 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); |
1069 | |
1070 | // skip errors for now |
1071 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
1072 | { |
1073 | chi2Result->SetBinError(j, 0); |
1074 | bayesResult->SetBinError(j, 0); |
1075 | } |
1076 | |
1077 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); |
1078 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
1079 | |
1080 | chi2Result->DrawCopy(""); |
1081 | bayesResult->DrawCopy("SAME"); |
1082 | } |
1083 | |
1084 | canvas3->cd(1); |
1085 | legend->Draw(); |
1086 | |
cfc19dd5 |
1087 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); |
dd701109 |
1088 | |
1089 | TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400); |
1090 | canvas2->Divide(2, 1); |
1091 | |
1092 | canvas2->cd(1); |
1093 | fitResultsChi2->SetMarkerStyle(20); |
1094 | fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); |
1095 | fitResultsChi2->Draw("AP"); |
1096 | |
1097 | fitResultsBayes->SetMarkerStyle(3); |
1098 | fitResultsBayes->SetMarkerColor(2); |
1099 | fitResultsBayes->Draw("P SAME"); |
1100 | |
1101 | gPad->SetLogy(); |
1102 | |
1103 | canvas2->cd(2); |
1104 | fitResultsChi2Limit->SetMarkerStyle(20); |
1105 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); |
1106 | fitResultsChi2Limit->Draw("AP"); |
1107 | |
1108 | fitResultsBayesLimit->SetMarkerStyle(3); |
1109 | fitResultsBayesLimit->SetMarkerColor(2); |
1110 | fitResultsBayesLimit->Draw("P SAME"); |
1111 | |
1112 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); |
1113 | canvas3->SaveAs(Form("%s.gif", canvas3->GetName())); |
cfc19dd5 |
1114 | } |
1115 | |
447c325d |
1116 | void DifferentSamples(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3) |
1117 | { |
1118 | gSystem->Load("libPWG0base"); |
1119 | |
1120 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1121 | |
1122 | TFile::Open(fileNameMC); |
1123 | mult->LoadHistograms("Multiplicity"); |
1124 | |
1125 | const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" }; |
1126 | |
1127 | TGraph* fitResultsChi2 = new TGraph; |
1128 | fitResultsChi2->SetTitle(";Input Dist ID;Chi2"); |
1129 | TGraph* fitResultsBayes = new TGraph; |
1130 | fitResultsBayes->SetTitle(";Input Dist ID;Chi2"); |
1131 | TGraph* fitResultsChi2Limit = new TGraph; |
1132 | fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach"); |
1133 | TGraph* fitResultsBayesLimit = new TGraph; |
1134 | fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach"); |
1135 | |
1136 | TCanvas* canvasA = new TCanvas("DifferentSamplesA", "DifferentSamplesA", 1200, 600); |
1137 | canvasA->Divide(4, 2); |
1138 | |
1139 | TCanvas* canvasB = new TCanvas("DifferentSamplesB", "DifferentSamplesB", 1200, 600); |
1140 | canvasB->Divide(4, 2); |
1141 | |
1142 | TCanvas* canvas4 = new TCanvas("DifferentSamples4", "DifferentSamples4", 1000, 400); |
1143 | canvas4->Divide(2, 1); |
1144 | |
1145 | TCanvas* canvas3 = new TCanvas("DifferentSamples3", "DifferentSamples3", 1000, 400); |
1146 | canvas3->Divide(2, 1); |
1147 | |
1148 | Float_t min = 1e10; |
1149 | Float_t max = 0; |
1150 | |
1151 | TH1* firstChi = 0; |
1152 | TH1* firstBayesian = 0; |
1153 | |
1154 | TLegend* legend = new TLegend(0.7, 0.7, 1, 1); |
1155 | |
1156 | TFile* file = TFile::Open("DifferentSamples.root", "RECREATE"); |
1157 | file->Close(); |
1158 | |
1159 | for (Int_t i=0; i<8; ++i) |
1160 | { |
1161 | TFile::Open(files[i]); |
1162 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2"); |
1163 | multESD->LoadHistograms("Multiplicity"); |
1164 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); |
1165 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i)); |
1166 | mc->Sumw2(); |
1167 | |
1168 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
1169 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); |
1170 | mult->DrawComparison(Form("DifferentSamples_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); |
1171 | |
1172 | Float_t chi2MC = 0; |
1173 | Int_t chi2MCLimit = 0; |
1174 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
1175 | fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC); |
1176 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit); |
1177 | min = TMath::Min(min, chi2MC); |
1178 | max = TMath::Max(max, chi2MC); |
1179 | |
1180 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i)); |
1181 | if (!firstChi) |
1182 | firstChi = (TH1*) chi2Result->Clone("firstChi"); |
1183 | |
1184 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100); |
1185 | mult->DrawComparison(Form("DifferentSamples_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
1186 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i)); |
1187 | if (!firstBayesian) |
1188 | firstBayesian = (TH1*) bayesResult->Clone("firstBayesian"); |
1189 | |
1190 | TFile* file = TFile::Open("DifferentSamples.root", "UPDATE"); |
1191 | mc->Write(); |
1192 | chi2Result->Write(); |
1193 | bayesResult->Write(); |
1194 | file->Close(); |
1195 | |
1196 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); |
1197 | fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC); |
1198 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit); |
1199 | |
1200 | min = TMath::Min(min, chi2MC); |
1201 | max = TMath::Max(max, chi2MC); |
1202 | mc->GetXaxis()->SetRangeUser(0, 150); |
1203 | chi2Result->GetXaxis()->SetRangeUser(0, 150); |
1204 | |
1205 | // skip errors for now |
1206 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
1207 | { |
1208 | chi2Result->SetBinError(j, 0); |
1209 | bayesResult->SetBinError(j, 0); |
1210 | } |
1211 | |
1212 | canvas4->cd(1); |
1213 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); |
1214 | tmp->SetTitle("Unfolded/MC;Npart;Ratio"); |
1215 | tmp->Divide(mc); |
1216 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1217 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1218 | tmp->SetLineColor(i+1); |
1219 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1220 | |
1221 | canvas4->cd(2); |
1222 | tmp = (TH1*) bayesResult->Clone("tmp"); |
1223 | tmp->SetTitle("Unfolded/MC;Npart;Ratio"); |
1224 | tmp->Divide(mc); |
1225 | tmp->SetLineColor(i+1); |
1226 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1227 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1228 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1229 | |
1230 | canvas3->cd(1); |
1231 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); |
1232 | tmp->SetTitle("Ratio to first result;Npart;Ratio"); |
1233 | tmp->Divide(firstChi); |
1234 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1235 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1236 | tmp->SetLineColor(i+1); |
1237 | legend->AddEntry(tmp, Form("%d", i)); |
1238 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1239 | |
1240 | canvas3->cd(2); |
1241 | tmp = (TH1*) bayesResult->Clone("tmp"); |
1242 | tmp->SetTitle("Ratio to first result;Npart;Ratio"); |
1243 | tmp->Divide(firstBayesian); |
1244 | tmp->SetLineColor(i+1); |
1245 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); |
1246 | tmp->GetXaxis()->SetRangeUser(0, 200); |
1247 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); |
1248 | |
1249 | if (i < 4) |
1250 | { |
1251 | canvasA->cd(i+1); |
1252 | } |
1253 | else |
1254 | canvasB->cd(i+1-4); |
1255 | |
1256 | mc->SetFillColor(kYellow); |
1257 | mc->DrawCopy(); |
1258 | chi2Result->SetLineColor(kRed); |
1259 | chi2Result->DrawCopy("SAME"); |
1260 | bayesResult->SetLineColor(kBlue); |
1261 | bayesResult->DrawCopy("SAME"); |
1262 | gPad->SetLogy(); |
1263 | |
1264 | if (i < 4) |
1265 | { |
1266 | canvasA->cd(i+5); |
1267 | } |
1268 | else |
1269 | canvasB->cd(i+5-4); |
1270 | |
1271 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); |
1272 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); |
1273 | |
1274 | // skip errors for now |
1275 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) |
1276 | { |
1277 | chi2Result->SetBinError(j, 0); |
1278 | bayesResult->SetBinError(j, 0); |
1279 | } |
1280 | |
1281 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); |
1282 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); |
1283 | |
1284 | chi2Result->DrawCopy(""); |
1285 | bayesResult->DrawCopy("SAME"); |
1286 | } |
1287 | |
1288 | canvas3->cd(1); |
1289 | legend->Draw(); |
1290 | |
1291 | canvasA->SaveAs(Form("%s.gif", canvasA->GetName())); |
1292 | canvasB->SaveAs(Form("%s.gif", canvasB->GetName())); |
1293 | |
1294 | TCanvas* canvas2 = new TCanvas("DifferentSamples2", "DifferentSamples2", 800, 400); |
1295 | canvas2->Divide(2, 1); |
1296 | |
1297 | canvas2->cd(1); |
1298 | fitResultsChi2->SetMarkerStyle(20); |
1299 | fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); |
1300 | fitResultsChi2->Draw("AP"); |
1301 | |
1302 | fitResultsBayes->SetMarkerStyle(3); |
1303 | fitResultsBayes->SetMarkerColor(2); |
1304 | fitResultsBayes->Draw("P SAME"); |
1305 | |
1306 | gPad->SetLogy(); |
1307 | |
1308 | canvas2->cd(2); |
1309 | fitResultsChi2Limit->SetMarkerStyle(20); |
1310 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); |
1311 | fitResultsChi2Limit->Draw("AP"); |
1312 | |
1313 | fitResultsBayesLimit->SetMarkerStyle(3); |
1314 | fitResultsBayesLimit->SetMarkerColor(2); |
1315 | fitResultsBayesLimit->Draw("P SAME"); |
1316 | |
1317 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); |
1318 | canvas3->SaveAs(Form("%s.gif", canvas3->GetName())); |
1319 | canvas4->SaveAs(Form("%s.gif", canvas4->GetName())); |
1320 | } |
1321 | |
cfc19dd5 |
1322 | void Merge(Int_t n, const char** files, const char* output) |
1323 | { |
447c325d |
1324 | // const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" }; |
1325 | |
1326 | |
cfc19dd5 |
1327 | gSystem->Load("libPWG0base"); |
1328 | |
1329 | AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n]; |
1330 | TList list; |
1331 | for (Int_t i=0; i<n; ++i) |
1332 | { |
1333 | TString name("Multiplicity"); |
1334 | if (i > 0) |
1335 | name.Form("Multiplicity%d", i); |
1336 | |
1337 | TFile::Open(files[i]); |
1338 | data[i] = new AliMultiplicityCorrection(name, name); |
1339 | data[i]->LoadHistograms("Multiplicity"); |
1340 | if (i > 0) |
1341 | list.Add(data[i]); |
1342 | } |
1343 | |
1344 | data[0]->Merge(&list); |
1345 | |
447c325d |
1346 | //data[0]->DrawHistograms(); |
cfc19dd5 |
1347 | |
1348 | TFile::Open(output, "RECREATE"); |
1349 | data[0]->SaveHistograms(); |
1350 | gFile->Close(); |
1351 | } |
1352 | |
1353 | void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root") |
1354 | { |
1355 | gSystem->Load("libPWG0base"); |
1356 | |
1357 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1358 | |
1359 | TFile::Open(fileName); |
1360 | mult->LoadHistograms("Multiplicity"); |
1361 | |
1362 | TF1* func = 0; |
1363 | |
1364 | if (caseNo >= 4) |
1365 | { |
447c325d |
1366 | 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, 250); |
cfc19dd5 |
1367 | func->SetParNames("scaling", "averagen", "k"); |
1368 | } |
1369 | |
1370 | switch (caseNo) |
1371 | { |
1372 | case 0: func = new TF1("flat", "1"); break; |
1373 | case 1: func = new TF1("flat", "501-x"); break; |
1374 | case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break; |
1375 | case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break; |
dd701109 |
1376 | case 4: func->SetParameters(1e7, 10, 2); break; |
447c325d |
1377 | case 5: func->SetParameters(1, 13, 7); break; |
dd701109 |
1378 | case 6: func->SetParameters(1e7, 30, 4); break; |
447c325d |
1379 | case 7: func->SetParameters(1e7, 30, 2); break; |
dd701109 |
1380 | case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break; |
cfc19dd5 |
1381 | |
1382 | default: return; |
1383 | } |
1384 | |
447c325d |
1385 | new TCanvas; |
1386 | func->Draw(); |
1387 | |
1388 | mult->SetGenMeasFromFunc(func, 3); |
cfc19dd5 |
1389 | |
dd701109 |
1390 | TFile::Open("out.root", "RECREATE"); |
1391 | mult->SaveHistograms(); |
1392 | |
1393 | //mult->ApplyBayesianMethod(2, kFALSE); |
1394 | //mult->ApplyMinuitFit(2, kFALSE); |
cfc19dd5 |
1395 | //mult->ApplyGaussianMethod(2, kFALSE); |
447c325d |
1396 | //mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx); |
dd701109 |
1397 | } |
1398 | |
1399 | void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2) |
1400 | { |
1401 | gSystem->Load("libPWG0base"); |
1402 | |
1403 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1404 | |
1405 | TFile::Open(fileName); |
1406 | mult->LoadHistograms("Multiplicity"); |
1407 | |
1408 | // empty under/overflow bins in x, otherwise Project3D takes them into account |
1409 | TH3* corr = mult->GetCorrelation(corrMatrix); |
1410 | for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j) |
1411 | { |
1412 | for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k) |
1413 | { |
1414 | corr->SetBinContent(0, j, k, 0); |
1415 | corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0); |
1416 | } |
1417 | } |
1418 | |
1419 | TH2* proj = (TH2*) corr->Project3D("zy"); |
1420 | |
1421 | // normalize correction for given nPart |
1422 | for (Int_t i=1; i<=proj->GetNbinsX(); ++i) |
1423 | { |
1424 | Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY()); |
1425 | if (sum <= 0) |
1426 | continue; |
1427 | |
1428 | for (Int_t j=1; j<=proj->GetNbinsY(); ++j) |
1429 | { |
1430 | // npart sum to 1 |
1431 | proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum); |
1432 | proj->SetBinError(i, j, proj->GetBinError(i, j) / sum); |
1433 | } |
1434 | } |
1435 | |
1436 | new TCanvas; |
447c325d |
1437 | proj->DrawCopy("COLZ"); |
dd701109 |
1438 | |
1439 | TH1* scaling = proj->ProjectionY("scaling", 1, 1); |
1440 | scaling->Reset(); |
1441 | scaling->SetMarkerStyle(3); |
1442 | //scaling->GetXaxis()->SetRangeUser(0, 50); |
1443 | TH1* mean = (TH1F*) scaling->Clone("mean"); |
1444 | TH1* width = (TH1F*) scaling->Clone("width"); |
1445 | |
1446 | TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500); |
1447 | lognormal->SetParNames("scaling", "mean", "sigma"); |
1448 | lognormal->SetParameters(1, 1, 1); |
1449 | lognormal->SetParLimits(0, 1, 1); |
1450 | lognormal->SetParLimits(1, 0, 100); |
1451 | lognormal->SetParLimits(2, 1e-3, 1); |
1452 | |
1453 | 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); |
1454 | nbd->SetParNames("scaling", "averagen", "k"); |
1455 | nbd->SetParameters(1, 13, 5); |
1456 | nbd->SetParLimits(0, 1, 1); |
1457 | nbd->SetParLimits(1, 1, 100); |
1458 | nbd->SetParLimits(2, 1, 1e8); |
1459 | |
1460 | TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50); |
1461 | poisson->SetParNames("scaling", "k", "deltax"); |
1462 | poisson->SetParameters(1, 1, 0); |
1463 | poisson->SetParLimits(0, 0, 10); |
1464 | poisson->SetParLimits(1, 0.01, 100); |
1465 | poisson->SetParLimits(2, 0, 10); |
1466 | |
1467 | TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50); |
1468 | mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth"); |
1469 | mygaus->SetParameters(1, 0, 1, 1, 0, 1); |
1470 | mygaus->SetParLimits(2, 1e-5, 10); |
1471 | mygaus->SetParLimits(4, 1, 1); |
1472 | mygaus->SetParLimits(5, 1e-5, 10); |
1473 | |
1474 | //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50); |
1475 | TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50); |
1476 | sqrt->SetParNames("ydelta", "exp", "xdelta"); |
1477 | sqrt->SetParameters(0, 0, 1); |
1478 | sqrt->SetParLimits(1, 0, 10); |
1479 | |
1480 | const char* fitWith = "gaus"; |
1481 | |
1482 | for (Int_t i=1; i<=150; ++i) |
1483 | { |
1484 | printf("Fitting %d...\n", i); |
1485 | |
1486 | TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e"); |
447c325d |
1487 | |
dd701109 |
1488 | //hist->GetXaxis()->SetRangeUser(0, 50); |
1489 | //lognormal->SetParameter(0, hist->GetMaximum()); |
1490 | hist->Fit(fitWith, "0 M", ""); |
1491 | |
1492 | TF1* func = hist->GetListOfFunctions()->FindObject(fitWith); |
1493 | |
447c325d |
1494 | if (0 && (i % 5 == 0)) |
dd701109 |
1495 | { |
447c325d |
1496 | pad = new TCanvas; |
dd701109 |
1497 | hist->Draw(); |
1498 | func->Clone()->Draw("SAME"); |
447c325d |
1499 | pad->SetLogy(); |
dd701109 |
1500 | } |
1501 | |
1502 | scaling->Fill(i, func->GetParameter(0)); |
1503 | mean->Fill(i, func->GetParameter(1)); |
1504 | width->Fill(i, func->GetParameter(2)); |
1505 | } |
1506 | |
1507 | TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500); |
1508 | log->SetParameters(0, 1, 1); |
1509 | log->SetParLimits(1, 0, 100); |
1510 | log->SetParLimits(2, 1e-3, 10); |
1511 | |
1512 | TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500); |
1513 | over->SetParameters(0, 1, 0); |
1514 | //over->SetParLimits(0, 0, 100); |
1515 | over->SetParLimits(1, 1e-3, 10); |
1516 | over->SetParLimits(2, 0, 100); |
1517 | |
1518 | c1 = new TCanvas("fitparams", "fitparams", 1200, 400); |
1519 | c1->Divide(3, 1); |
1520 | |
1521 | c1->cd(1); |
1522 | scaling->Draw("P"); |
1523 | |
1524 | //TF1* scalingFit = new TF1("mypol0", "[0]"); |
1525 | TF1* scalingFit = over; |
447c325d |
1526 | scaling->Fit(scalingFit, "", "", 3, 140); |
1527 | scalingFit->SetRange(0, 200); |
1528 | scalingFit->Draw("SAME"); |
dd701109 |
1529 | |
1530 | c1->cd(2); |
1531 | mean->Draw("P"); |
1532 | |
1533 | //TF1* meanFit = log; |
1534 | TF1* meanFit = new TF1("mypol1", "[0]+[1]*x"); |
447c325d |
1535 | mean->Fit(meanFit, "", "", 3, 140); |
1536 | meanFit->SetRange(0, 200); |
1537 | meanFit->Draw("SAME"); |
dd701109 |
1538 | |
1539 | c1->cd(3); |
1540 | width->Draw("P"); |
1541 | |
1542 | //TF1* widthFit = over; |
447c325d |
1543 | TF1* widthFit = new TF1("mypol", "[0]+[1]*TMath::Sqrt([2]*x)"); |
1544 | widthFit->SetParLimits(2, 1e-5, 1e5); |
1545 | width->Fit(widthFit, "", "", 5, 140); |
1546 | widthFit->SetRange(0, 200); |
1547 | widthFit->Draw("SAME"); |
dd701109 |
1548 | |
1549 | // build new correction matrix |
1550 | TH2* new = (TH2*) proj->Clone("new"); |
1551 | new->Reset(); |
1552 | Float_t x, y; |
1553 | for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) |
1554 | { |
1555 | TF1* func = (TF1*) gROOT->FindObject(fitWith); |
1556 | x = new->GetXaxis()->GetBinCenter(i); |
1557 | //if (i == 1) |
1558 | // x = 0.1; |
1559 | x++; |
1560 | func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); |
1561 | printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); |
1562 | |
1563 | for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) |
1564 | { |
447c325d |
1565 | if (i < 11) |
dd701109 |
1566 | { |
1567 | // leave bins 1..20 untouched |
1568 | new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j)); |
1569 | } |
1570 | else |
1571 | { |
1572 | y = new->GetYaxis()->GetBinCenter(j); |
1573 | if (j == 1) |
1574 | y = 0.1; |
1575 | if (func->Eval(y) > 1e-4) |
1576 | new->SetBinContent(i, j, func->Eval(y)); |
1577 | } |
1578 | } |
1579 | } |
1580 | |
1581 | // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0 |
1582 | // we take the values from the old response matrix |
1583 | //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) |
1584 | // new->SetBinContent(i, 1, proj->GetBinContent(i, 1)); |
1585 | |
1586 | //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) |
1587 | // new->SetBinContent(1, j, proj->GetBinContent(1, j)); |
1588 | |
1589 | // normalize correction for given nPart |
1590 | for (Int_t i=1; i<=new->GetNbinsX(); ++i) |
1591 | { |
1592 | Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY()); |
1593 | if (sum <= 0) |
1594 | continue; |
1595 | |
1596 | for (Int_t j=1; j<=new->GetNbinsY(); ++j) |
1597 | { |
1598 | // npart sum to 1 |
1599 | new->SetBinContent(i, j, new->GetBinContent(i, j) / sum); |
1600 | new->SetBinError(i, j, new->GetBinError(i, j) / sum); |
1601 | } |
1602 | } |
1603 | |
1604 | new TCanvas; |
1605 | new->Draw("COLZ"); |
1606 | |
1607 | TH2* diff = (TH2*) new->Clone("diff"); |
1608 | diff->Add(proj, -1); |
1609 | |
1610 | new TCanvas; |
1611 | diff->Draw("COLZ"); |
1612 | diff->SetMinimum(-0.05); |
1613 | diff->SetMaximum(0.05); |
1614 | |
1615 | corr->Reset(); |
1616 | |
1617 | for (Int_t i=1; i<=new->GetNbinsX(); ++i) |
1618 | for (Int_t j=1; j<=new->GetNbinsY(); ++j) |
1619 | corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2, i, j, new->GetBinContent(i, j)); |
1620 | |
1621 | new TCanvas; |
1622 | corr->Project3D("zy")->Draw("COLZ"); |
1623 | |
1624 | TFile::Open("out.root", "RECREATE"); |
1625 | mult->SaveHistograms(); |
1626 | |
1627 | TH1* proj1 = proj->ProjectionY("proj1", 36, 36); |
1628 | TH1* proj2 = new->ProjectionY("proj2", 36, 36); |
1629 | proj2->SetLineColor(2); |
1630 | |
1631 | new TCanvas; |
1632 | proj1->Draw(); |
1633 | proj2->Draw("SAME"); |
cfc19dd5 |
1634 | } |
447c325d |
1635 | |
1636 | void GetCrossSections(const char* fileName) |
1637 | { |
1638 | gSystem->Load("libPWG0base"); |
1639 | |
1640 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); |
1641 | |
1642 | TFile::Open(fileName); |
1643 | mult->LoadHistograms("Multiplicity"); |
1644 | |
1645 | TH1* xSection2 = mult->GetCorrelation(3)->Project3D("y")->Clone("xSection2"); |
1646 | xSection2->Sumw2(); |
1647 | xSection2->Scale(1.0 / xSection2->Integral()); |
1648 | |
1649 | TH1* xSection15 = mult->GetCorrelation(2)->Project3D("y")->Clone("xSection15"); |
1650 | xSection15->Sumw2(); |
1651 | xSection15->Scale(1.0 / xSection15->Integral()); |
1652 | |
1653 | TFile::Open("crosssection.root", "RECREATE"); |
1654 | xSection2->Write(); |
1655 | xSection15->Write(); |
1656 | gFile->Close(); |
1657 | } |