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