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