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