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