25fd74fee7a8f30dee7e120049733f53d8ee3e01
[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     AliUnfolding::SetChi2Regularization(AliUnfolding::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     //mult->ApplyNBDFit(histID, fullPhaseSpace, eventType);
103   }
104   else
105   {
106     mult->ApplyBayesianMethod(histID, fullPhaseSpace, eventType, 1, 10, 0, kTRUE);
107   }
108
109   mult->SetMultiplicityMC(histID, eventType, hist2);
110   
111   Printf("Writing result to %s", targetFile);
112   TFile* file = TFile::Open(targetFile, "RECREATE");
113   mult->SaveHistograms();
114   file->Write();
115   file->Close();
116
117   mult->DrawComparison((chi2) ? "MinuitChi2" : "Bayesian", histID, fullPhaseSpace, kTRUE, mcCompare);
118 }
119
120 void 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");
236 }
237
238 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)
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
252   histRAW = chi2->GetMultiplicityESD(histID)->ProjectionY("raw", 1, chi2->GetMultiplicityESD(histID)->GetNbinsX());
253   histRAW->Scale(1.0 / histRAW->Integral());
254
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();
262   c->SetGridx();
263   c->SetGridy();
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);
271   histRAW->SetLineColor(3);
272
273   histC->DrawCopy("HISTE");
274   histB->DrawCopy("HISTE SAME");
275   histRAW->DrawCopy("SAME");
276
277   legend = new TLegend(0.2, 0.2, 0.4, 0.4);
278   legend->SetFillColor(0);
279
280   legend->AddEntry(histC, label1);
281   legend->AddEntry(histB, label2);
282   legend->AddEntry(histRAW, "raw ESD");
283
284   histGraph = 0;
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)
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());
327     histMC->Sumw2();
328     histMC->Scale(1.0 / histMC->Integral());
329
330     histMC->Draw("HISTE SAME");
331     histMC->SetLineColor(4);
332     legend->AddEntry(histMC, "MC");
333   }
334
335   legend->Draw();
336
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));
357
358       /*
359       if (histGraph)
360       {
361         histGraph->SetBinContent(bin, histGraph->GetBinContent(bin) / histMC->GetBinContent(bin));
362         histGraph->SetBinError(bin, 0);
363       }
364       */
365
366       // TODO errors?
367       histC->SetBinError(bin, 0);
368       histB->SetBinError(bin, 0);
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");
377   /*
378   if (histGraph)
379     histGraph->Draw("HIST SAME");
380   */
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
404 void 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
445   c->SaveAs(Form("%s.gif", c->GetName()));
446   c->SaveAs(Form("%s.eps", c->GetName()));
447 }
448
449 void* 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"));
487
488   return mult;
489 }
490
491 const 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;
500     case AliMultiplicityCorrection::kLog   :    return "Log"; break;
501   }
502   return 0;
503 }
504
505 const 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;
512   }
513   return 0;
514 }
515
516 void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", const char* targetDir = "bayesian", Int_t histID = 1)
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
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};
537
538   for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
539   //for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
540   {
541     TString tmp;
542     tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
543
544     TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600);
545
546     for (Int_t i = 1; i <= 5; i++)
547     {
548       Int_t iterArray[5] = {2, 5, 10, 20, -1};
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;
555         fitResultsMC[region]->SetTitle(Form("%d iterations - reg. %d", iter, region+1));
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);
560         //fitResultsMC[region]->SetMarkerStyle(markers[(i-1) * AliMultiplicityCorrection::kQualityRegions + region]);
561         fitResultsMC[region]->SetMarkerStyle(markers[i-1]);
562         fitResultsMC[region]->SetLineColor(colors[i-1]);
563         fitResultsMC[region]->SetMarkerColor(colors[i-1]);
564       }
565
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
572       fitResultsRes->SetFillColor(0);
573       fitResultsRes->SetMarkerStyle(markers[i-1]);
574       fitResultsRes->SetMarkerColor(colors[i-1]);
575       fitResultsRes->SetLineColor(colors[i-1]);
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
582         mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter, 0, kFALSE);
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
586         for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
587           fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
588
589         fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
590       }
591
592       graphFile->cd();
593       for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
594         fitResultsMC[region]->Write();
595
596       fitResultsRes->Write();
597     }
598   }
599
600   graphFile->Close();
601 }
602
603 void 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)
622   {
623     canvas->SetLogx();
624     canvas->SetLogy();
625   }
626   canvas->SetTopMargin(0.05);
627   canvas->SetGridx();
628   canvas->SetGridy();
629
630   TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98);
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
641   Float_t yMinRegion[3];
642   for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
643     yMinRegion[i] = 1e20;
644
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
652     if (!mc)
653       break;
654
655     xaxis = mc->GetXaxis()->GetTitle();
656     yaxis = mc->GetYaxis()->GetTitle();
657
658     mc->Print();
659
660     if (res)
661       res->Print();
662
663     xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
664     yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
665
666     xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
667     yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
668
669     if (plotRes && res)
670     {
671       xMin = TMath::Min(xMin, res->GetXaxis()->GetXmin());
672       yMin = TMath::Min(yMin, res->GetYaxis()->GetXmin());
673
674       xMax = TMath::Max(xMax, res->GetXaxis()->GetXmax());
675       yMax = TMath::Max(yMax, res->GetYaxis()->GetXmax());
676     }
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++;
682   }
683
684   for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
685     Printf("Minimum for region %d is %f", i, yMinRegion[i]);
686
687   if (type >= 0)
688   {
689     xaxis = "smoothing parameter";
690   }
691   else if (type == -1)
692   {
693     xaxis = "weight parameter";
694     xMax *= 5;
695   }
696   //yaxis = "P_{1} (2 <= t <= 150)";
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");
707   dummy->GetYaxis()->SetMoreLogLabels(1);
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
716     //printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res);
717
718     if (!mc)
719       break;
720
721     printf("Loaded %d sucessful.\n", count);
722
723     if (type == -1)
724     {
725       legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3));
726     }
727     else
728       legend->AddEntry(mc);
729
730     mc->Draw("SAME PC");
731
732     if (plotRes && res)
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
747 void EvaluateDrawResultRegions(const char* targetDir, Int_t type = 0)
748 {
749   loadlibs();
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
764   TCanvas* canvas = new TCanvas(name, name, 1200, 800);
765   //canvas->Divide(1, AliMultiplicityCorrection::kQualityRegions, 0, 0);
766   canvas->SetTopMargin(0);
767   canvas->SetBottomMargin(0);
768   canvas->SetRightMargin(0.05);
769   canvas->Range(0, 0, 1, 1);
770
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); 
776
777   Float_t yMinRegion[3];
778   for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
779     yMinRegion[i] = 1e20;
780
781   for (Int_t region = 0; region <= AliMultiplicityCorrection::kQualityRegions; region++)
782   {
783     canvas->cd();
784     pad[region]->Draw();
785     pad[region]->cd();
786     pad[region]->SetRightMargin(0.05);
787     
788     if (type == -1)
789     {
790       pad[region]->SetLogx();
791     }
792     pad[region]->SetLogy();
793     
794     pad[region]->SetGridx();
795     pad[region]->SetGridy();
796
797     TLegend* legend = new TLegend(0.5, 0.4, 0.85, 0.85);
798     legend->SetFillColor(0);
799     //legend->SetNColumns(3);
800     legend->SetTextSize(0.06);
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
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       }
833
834       xaxis = mc->GetXaxis()->GetTitle();
835       yaxis = mc->GetYaxis()->GetTitle();
836
837       //mc->Print();
838
839       xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
840       yMin = TMath::Min(yMin, TMath::MinElement(mc->GetN(), mc->GetY()));
841
842       xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
843       yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
844       
845       //Printf("%f %f %f %f", xMin, xMax, yMin, yMax);
846     }
847
848     if (type >= 0)
849     {
850       xaxis = "Smoothing parameter #alpha";
851     }
852     else if (type == -1)
853     {
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)";
860     }
861     else
862       yaxis = "Q_{2} ";
863
864     printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
865
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);
884     dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
885
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    
896
897     count = 0;
898
899     while (1)
900     {
901       count++;
902
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
912       if (!mc)
913         break;
914       //printf("Loaded %d sucessful.\n", count);
915
916       if (type == -1)
917       {
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");
922       }
923       else
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       }
937
938       mc->Draw("SAME PC");
939
940     }
941
942     if (region == 2)
943       legend->Draw();
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
954 void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", const char* targetDir = "chi2compare", Int_t histID = 1)
955 {
956   gSystem->mkdir(targetDir);
957
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
971   Int_t count = 0; // just to order the saved images...
972   Int_t colors[] = {1, 2, 4, 3};
973   Int_t markers[] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3};
974
975   TGraph* fitResultsRes = 0;
976
977   TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE");
978
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)
982   {
983     TGraph* fitResultsMC[3];
984     for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
985     {
986       fitResultsMC[region] = new TGraph;
987       fitResultsMC[region]->SetTitle(Form("Eq. (%d) - reg. %d", type+9, region+1));
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);
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]);
997     }
998
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");
1003
1004     fitResultsRes->SetFillColor(0);
1005     fitResultsRes->SetMarkerStyle(markers[type-1]);
1006     fitResultsRes->SetMarkerColor(colors[type-1]);
1007     fitResultsRes->SetLineColor(colors[type-1]);
1008
1009     for (Int_t i=0; i<15; ++i)
1010     {
1011       Float_t weight = TMath::Power(TMath::Sqrt(10), i+1);
1012       //Float_t weight = TMath::Power(10, i+2);
1013
1014       //if (type == AliMultiplicityCorrection::kEntropy) weight = 1e4 * (i+1) * 1.5;
1015
1016       Float_t chi2MC = 0;
1017       Float_t residuals = 0;
1018       Float_t chi2Limit = 0;
1019
1020       TString runName;
1021       runName.Form("MinuitChi2_%02d_%d_%f", count++, type, weight);
1022
1023       mult->SetRegularizationParameters(type, weight);
1024       Int_t status = mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
1025       mult->DrawComparison(Form("%s/%s", targetDir, runName.Data()), histID, kFALSE, kTRUE, hist2->ProjectionY("mymc", 1, hist2->GetNbinsX()));
1026       if (status != 0)
1027       {
1028         printf("MINUIT did not succeed. Skipping...\n");
1029         continue;
1030       }
1031
1032       mult->GetComparisonResults(&chi2MC, 0, &residuals);
1033       TH1* result = mult->GetMultiplicityESDCorrected(histID);
1034       result->SetName(runName);
1035       result->Write();
1036
1037       for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1038         fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
1039
1040       fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
1041     }
1042
1043     graphFile->cd();
1044     for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1045       fitResultsMC[region]->Write();
1046     fitResultsRes->Write();
1047   }
1048
1049   graphFile->Close();
1050 }
1051
1052 void 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
1061 void 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
1070 void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
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
1087   TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 1200);
1088   canvas->Divide(3, 3);
1089
1090   Int_t count = 0;
1091
1092   for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
1093   {
1094     TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY("mymc");
1095     mc->Sumw2();
1096
1097     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
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
1109 /*    // skip errors for now
1110     for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
1111     {
1112       chi2Result->SetBinError(i, 0);
1113       bayesResult->SetBinError(i, 0);
1114     }*/
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);
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, "");
1130
1131     chi2ResultRatio->GetYaxis()->SetRangeUser(0.5, 1.5);
1132
1133     chi2ResultRatio->DrawCopy("HIST");
1134     bayesResultRatio->DrawCopy("SAME HIST");
1135
1136     canvas->cd(count + 6);
1137     chi2Result->Divide(chi2Result, bayesResult, 1, 1, "");
1138     chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1139     chi2Result->DrawCopy("HIST");
1140   }
1141
1142   canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
1143   canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
1144 }
1145
1146 void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
1147 {
1148   gSystem->Load("libPWG0base");
1149
1150   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1151
1152   TFile::Open(fileNameMC);
1153   mult->LoadHistograms("Multiplicity");
1154
1155   const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" };
1156
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
1174   TGraph* fitResultsChi2Limit = new TGraph;  fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
1175   TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
1176   TGraph* fitResultsChi2Res = new TGraph;       fitResultsChi2Res->SetTitle(";Nevents;Chi2");
1177   TGraph* fitResultsBayesRes = new TGraph;      fitResultsBayesRes->SetTitle(";Nevents;Chi2");
1178
1179   fitResultsChi2Limit->SetName("fitResultsChi2Limit");
1180   fitResultsBayesLimit->SetName("fitResultsBayesLimit");
1181   fitResultsChi2Res->SetName("fitResultsChi2Res");
1182   fitResultsBayesRes->SetName("fitResultsBayesRes");
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
1190   TFile* file = TFile::Open("StatisticsPlot.root", "RECREATE");
1191   file->Close();
1192
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();
1201     TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i));
1202
1203     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1204     mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
1205     mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
1206
1207     Int_t chi2MCLimit = 0;
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     }
1216     fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit);
1217     fitResultsChi2Res->SetPoint(fitResultsChi2Res->GetN(), nEntries, chi2Residuals);
1218
1219     TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
1220
1221     mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0, kFALSE);
1222     mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
1223     TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
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     }
1231     fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit);
1232     fitResultsBayesRes->SetPoint(fitResultsBayesRes->GetN(), nEntries, chi2Residuals);
1233
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");
1269
1270     TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
1271     mc->Write();
1272     chi2Result->Write();
1273     bayesResult->Write();
1274     file->Close();
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);
1284
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   }
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()));
1306
1307   TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
1308
1309   for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1310   {
1311     fitResultsChi2[region]->Write();
1312     fitResultsBayes[region]->Write();
1313   }
1314   fitResultsChi2Limit->Write();
1315   fitResultsBayesLimit->Write();
1316   fitResultsChi2Res->Write();
1317   fitResultsBayesRes->Write();
1318   file->Close();
1319 }
1320
1321 void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 1)
1322 {
1323   loadlibs();
1324
1325   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1326
1327   TFile::Open(fileNameMC);
1328   mult->LoadHistograms("Multiplicity");
1329
1330   const char* files[] = { "multiplicityESD.root", "multiplicityESD_100k_1.root", "multiplicityESD_100k_2.root", "multiplicityESD_100k_1.root", "multiplicityESD_100k_2.root" }
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));
1337   TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY("mc", 1, 1);
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
1363   TFile* file = TFile::Open("StartingConditions.root", "RECREATE");
1364   mc->Write();
1365   file->Close();
1366
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     {
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       
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
1399     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 1e5);
1400     mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond);
1401     //mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
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
1411     TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
1412     if (!firstChi)
1413       firstChi = (TH1*) chi2Result->Clone("firstChi");
1414
1415     mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 10, startCond);
1416     //mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
1417     TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
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
1425     TFile* file = TFile::Open("StartingConditions.root", "UPDATE");
1426     chi2Result->Write();
1427     bayesResult->Write();
1428     file->Close();
1429
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
1491   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
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()));
1518 }
1519
1520 void Merge(Int_t n, const char** files, const char* output)
1521 {
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
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
1544   //data[0]->DrawHistograms();
1545
1546   TFile::Open(output, "RECREATE");
1547   data[0]->SaveHistograms();
1548   gFile->Close();
1549 }
1550
1551 void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
1552 {
1553   loadlibs();
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   {
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);
1566     func->SetParNames("scaling", "averagen", "k");
1567   }
1568
1569   switch (caseNo)
1570   {
1571     case 0: func = new TF1("flat", "1000"); break;
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;
1575     case 4: func->SetParameters(2e5, 15, 2); break;
1576     case 5: func->SetParameters(1, 13, 7); break;
1577     case 6: func->SetParameters(1e7, 30, 4); break;
1578     case 7: func->SetParameters(1e7, 30, 2); break; // ***
1579     case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break;
1580
1581     default: return;
1582   }
1583
1584   new TCanvas;
1585   func->Draw();
1586
1587   mult->SetGenMeasFromFunc(func, 1);
1588
1589   TFile::Open("out.root", "RECREATE");
1590   mult->SaveHistograms();
1591
1592   new TCanvas; mult->GetMultiplicityESD(1)->ProjectionY()->DrawCopy();
1593   new TCanvas; mult->GetMultiplicityVtx(1)->ProjectionY("mc", 1, mult->GetMultiplicityVtx(1)->GetNbinsX())->DrawCopy();
1594
1595   //mult->ApplyBayesianMethod(2, kFALSE);
1596   //mult->ApplyMinuitFit(2, kFALSE);
1597   //mult->ApplyGaussianMethod(2, kFALSE);
1598   //mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
1599 }
1600
1601 void 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;
1639   proj->DrawCopy("COLZ");
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");
1689
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
1696     if (0 && (i % 5 == 0))
1697     {
1698       pad = new TCanvas;
1699       hist->Draw();
1700       func->Clone()->Draw("SAME");
1701       pad->SetLogy();
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;
1728   scaling->Fit(scalingFit, "", "", 3, 140);
1729   scalingFit->SetRange(0, 200);
1730   scalingFit->Draw("SAME");
1731
1732   c1->cd(2);
1733   mean->Draw("P");
1734
1735   //TF1* meanFit = log;
1736   TF1* meanFit = new TF1("mypol1", "[0]+[1]*x");
1737   mean->Fit(meanFit, "", "", 3, 140);
1738   meanFit->SetRange(0, 200);
1739   meanFit->Draw("SAME");
1740
1741   c1->cd(3);
1742   width->Draw("P");
1743
1744   //TF1* widthFit = over;
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");
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     {
1767       if (i < 11)
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");
1836 }
1837
1838 void 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
1886 void 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 }
1908
1909 void 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
1939 void BuildResponseFromTree(const char* fileName, const char* target)
1940 {
1941   //
1942   // builds several response matrices with different particle ratios (systematic study)
1943   // 
1944   // WARNING doesn't work uncompiled, see test.C
1945
1946   loadlibs();
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;
1956   Int_t secondaries = 0;
1957   Int_t doubleCount = 0;
1958
1959   for (Int_t num = 0; num < 9; num++)
1960   {
1961     TString name;
1962     name.Form("Multiplicity_%d", num);
1963     AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(name, name);
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;
1972       case 2 : ratio[1] = 1.5; break;
1973       case 3 : ratio[2] = 0.5; break;
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;
1977       case 7 : ratio[1] = 0.5; ratio[2] = 1.5; break;
1978       case 8 : ratio[1] = 1.5; ratio[2] = 0.5; break;
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
2001       // secondaries are already part of meas!
2002       secondaries += f[10];
2003
2004       // double counted are already part of meas!
2005       doubleCount += f[11];
2006
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!
2008       meas += f[9];
2009
2010       //Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]);
2011
2012       fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, 0, meas, meas, meas, meas);
2013       // HACK all as kND = 1
2014       fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, 1, gene, gene, gene, gene, 0);
2015       fMultiplicity->FillMeasured(f[0], meas, meas, meas, meas);
2016     }
2017
2018     //fMultiplicity->DrawHistograms();
2019
2020     TFile* file = TFile::Open(target, "UPDATE");
2021     fMultiplicity->SaveHistograms();
2022     file->Close();
2023
2024     if (num == 0)
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     }
2030   }
2031 }
2032
2033 void ParticleRatioStudy()
2034 {
2035   loadlibs();
2036
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
2045 void MergeModifyCrossSection(const char* output = "multiplicityMC_xsection.root")
2046 {
2047   const char* files[] = { "multiplicityMC_nd.root", "multiplicityMC_sd.root", "multiplicityMC_dd.root" };
2048
2049   loadlibs();
2050
2051   TFile::Open(output, "RECREATE");
2052   gFile->Close();
2053
2054   for (Int_t num = 0; num < 7; num++)
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]);
2089         data[i]->GetMultiplicityNSD(j)->Scale(ratio[i]);
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 }
2118
2119 void 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
2188 void 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 }