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