Fix for CMake
[u/mrichter/AliRoot.git] / PWG0 / multiplicity / plots.C
1 /* $Id$ */
2
3 //
4 // plots for the note about multplicity measurements
5 //
6
7 #if !defined(__CINT__) || defined(__MAKECINT__)
8
9 #include <TCanvas.h>
10 #include <TPad.h>
11 #include <TH1F.h>
12 #include <TH2F.h>
13 #include <TH3F.h>
14 #include <TLine.h>
15 #include <TF1.h>
16 #include <TSystem.h>
17 #include <TFile.h>
18 #include <TLegend.h>
19 #include <TStopwatch.h>
20 #include <TROOT.h>
21 #include <TGraph.h>
22 #include <TMath.h>
23 #include <TPaveText.h>
24 #include <TImage.h>
25 #include <TLatex.h>
26
27 #include "AliMultiplicityCorrection.h"
28 #include "AliCorrection.h"
29 #include "AliCorrectionMatrix3D.h"
30
31 #endif
32
33 const char* correctionFile = "multiplicityMC.root";
34 const char* measuredFile   = "multiplicityESD.root";
35 Int_t etaRange = 1;
36 Int_t displayRange = 80; // axis range
37 Int_t ratioRange = 151;   // range to calculate difference
38 Int_t longDisplayRange = 120;
39
40 const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
41 const char* measuredFileTPC   = "multiplicityMC_TPC_0.6M.root";
42 Int_t etaRangeTPC = 1;
43
44 void loadlibs()
45 {
46   gSystem->Load("libANALYSIS");
47   gSystem->Load("libPWG0base");
48 }
49
50 void SetTPC()
51 {
52   correctionFile = correctionFileTPC;
53   measuredFile = measuredFileTPC;
54   etaRange = etaRangeTPC;
55   displayRange = 100;
56   ratioRange = 76;
57   longDisplayRange = 100;
58 }
59
60 const char* GetMultLabel(Int_t etaR = -1, Bool_t trueM = kTRUE)
61 {
62         if (etaR == -1)
63                 etaR = etaRange;
64                 
65         TString tmpStr((trueM) ? "True " : "Measured ");
66
67         tmpStr += "multiplicity";
68         //return Form("%s", tmpStr.Data());
69                 
70         if (etaR == 4)
71         {
72                 tmpStr += " (full phase space)";
73         }
74         else
75                 tmpStr += Form(" in |#eta| < %.1f", (etaR+1)* 0.5);
76         return Form("%s", tmpStr.Data());
77 }
78
79 void Smooth(TH1* hist, Int_t windowWidth = 20)
80 {
81   TH1* clone = (TH1*) hist->Clone("clone");
82   for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
83   {
84     Int_t min = TMath::Max(2, bin-windowWidth);
85     Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
86     Float_t average = clone->Integral(min, max) / (max - min + 1);
87
88     hist->SetBinContent(bin, average);
89     hist->SetBinError(bin, 0);
90   }
91
92   delete clone;
93 }
94
95 void responseMatrixPlot(const char* fileName = 0)
96 {
97   loadlibs();
98
99   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
100
101   if (fileName == 0)
102     fileName = correctionFile;
103   
104   TFile::Open(fileName);
105   mult->LoadHistograms("Multiplicity");
106
107   // empty under/overflow bins in x, otherwise Project3D takes them into account
108   TH1* hist = mult->GetCorrelation(etaRange);
109   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
110   {
111     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
112     {
113       hist->SetBinContent(0, y, z, 0);
114       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
115     }
116   }
117   hist = ((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy");
118   hist->SetStats(kFALSE);
119
120   hist->SetTitle(Form(";%s;%s;Entries", GetMultLabel(), GetMultLabel(etaRange, kFALSE)));
121   hist->GetXaxis()->SetRangeUser(0, longDisplayRange);
122   hist->GetYaxis()->SetRangeUser(0, longDisplayRange);
123   
124   hist->GetYaxis()->SetTitleOffset(1.3);
125   hist->GetZaxis()->SetTitleOffset(1.2);
126
127   TCanvas* canvas = new TCanvas("c1", "c1", 600, 600);
128   canvas->SetRightMargin(0.15);
129   canvas->SetTopMargin(0.05);
130
131   gPad->SetLogz();
132   hist->Draw("COLZ");
133
134   canvas->SaveAs("responsematrix.eps");
135 }
136
137 void multPythiaPhojet()
138 {
139   loadlibs();
140   
141   TFile::Open("LHC08c11_10TeV_0.5T/mb1/spd/multiplicity.root");
142   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
143   mult->LoadHistograms("Multiplicity");
144   hist1 = mult->GetMultiplicityINEL(1)->ProjectionY();
145   hist1->Sumw2();
146   
147   TFile::Open("LHC08c15_10TeV_0.5T_Phojet/mb1/spd/multiplicity.root");
148   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
149   mult->LoadHistograms("Multiplicity");
150   hist2 = mult->GetMultiplicityINEL(1)->ProjectionY();
151   hist2->Sumw2();
152   
153   legend = new TLegend(0.6, 0.7, 0.9, 0.9);
154   legend->SetFillColor(0);
155   legend->SetTextSize(0.04);
156   legend->AddEntry(hist1, "Pythia", "L");
157   legend->AddEntry(hist2, "Phojet", "L");
158
159   c1 = new TCanvas("c", "c", 600, 600);
160   c1->SetTopMargin(0.05);
161   c1->SetRightMargin(0.05);
162   c1->SetLeftMargin(0.12);
163   c1->SetGridx();
164   c1->SetGridy(); 
165   c1->SetLogy();
166   
167   //hist1->SetMarkerStyle(20);
168   //hist2->SetMarkerStyle(24);
169   //hist2->SetMarkerColor(2);
170   hist1->SetLineWidth(2);
171   hist2->SetLineWidth(2);
172   hist2->SetLineStyle(2);
173   hist2->SetLineColor(2);
174   
175   hist1->Scale(1.0 / hist1->Integral());
176   hist2->Scale(1.0 / hist2->Integral());
177   
178   hist1->SetStats(0);
179   hist1->GetYaxis()->SetTitleOffset(1.3);
180   hist1->GetXaxis()->SetRangeUser(0, 100);
181   hist1->SetTitle(";N_{ch};P(N_{ch})");
182   
183   hist1->Draw("");
184   hist2->Draw("SAME");
185   legend->Draw();
186   
187   c1->SaveAs("mult_pythia_phojet.eps");
188 }
189
190 TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
191 {
192   // normalize unfolded result to mc hist
193   result->Scale(1.0 / result->Integral(2, displayRange));
194   result->Scale(mcHist->Integral(2, displayRange));
195
196   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
197   canvas->Range(0, 0, 1, 1);
198
199   TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
200   pad1->Draw();
201
202   TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
203   pad2->Draw();
204
205   pad1->SetRightMargin(0.05);
206   pad2->SetRightMargin(0.05);
207
208   // no border between them
209   pad1->SetBottomMargin(0);
210   pad2->SetTopMargin(0);
211
212   pad1->cd();
213   pad1->SetGridx();
214   pad1->SetGridy();
215
216   mcHist->GetXaxis()->SetLabelSize(0.06);
217   mcHist->GetYaxis()->SetLabelSize(0.06);
218   mcHist->GetXaxis()->SetTitleSize(0.06);
219   mcHist->GetYaxis()->SetTitleSize(0.06);
220   mcHist->GetYaxis()->SetTitleOffset(0.6);
221
222   mcHist->GetXaxis()->SetRangeUser(0, displayRange);
223
224   mcHist->SetTitle(Form(";%s;Entries", GetMultLabel()));
225   mcHist->SetStats(kFALSE);
226
227   mcHist->DrawCopy("HIST E");
228   gPad->SetLogy();
229
230   result->SetLineColor(2);
231   result->SetMarkerColor(2);
232   result->SetMarkerStyle(5);
233   result->DrawCopy("SAME PE");
234
235   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
236   legend->AddEntry(mcHist, "True distribution");
237   legend->AddEntry(result, "Unfolded distribution", "P");
238   legend->SetFillColor(0);
239   legend->Draw();
240
241   pad2->cd();
242   pad2->SetBottomMargin(0.15);
243
244   // calculate ratio
245   mcHist->Sumw2();
246   TH1* ratio = (TH1*) mcHist->Clone("ratio");
247   result->Sumw2();
248   ratio->Divide(ratio, result, 1, 1, "");
249   ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
250   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
251
252   ratio->DrawCopy();
253
254   // get average of ratio
255   Float_t sum = 0;
256   for (Int_t i=2; i<=ratioRange; ++i)
257   {
258     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
259   }
260   sum /= ratioRange-1;
261
262   printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
263
264   TLine* line = new TLine(0, 1, displayRange, 1);
265   line->SetLineWidth(2);
266   line->Draw();
267
268   line = new TLine(0, 1.1, displayRange, 1.1);
269   line->SetLineWidth(2);
270   line->SetLineStyle(2);
271   line->Draw();
272   line = new TLine(0, 0.9, displayRange, 0.9);
273   line->SetLineWidth(2);
274   line->SetLineStyle(2);
275   line->Draw();
276
277   canvas->Modified();
278
279   canvas->SaveAs(epsName);
280
281   return canvas;
282 }
283
284 TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
285 {
286   // draws the 3 plots in the upper plot
287   // draws the ratio between result1 and result2 in the lower plot
288
289   // normalize unfolded result to mc hist
290   result1->Scale(1.0 / result1->Integral(2, displayRange));
291   result1->Scale(mcHist->Integral(2, displayRange));
292   result2->Scale(1.0 / result2->Integral(2, displayRange));
293   result2->Scale(mcHist->Integral(2, displayRange));
294
295   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
296   canvas->Range(0, 0, 1, 1);
297
298   TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
299   pad1->Draw();
300
301   TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
302   pad2->Draw();
303
304   pad1->SetRightMargin(0.05);
305   pad2->SetRightMargin(0.05);
306
307   // no border between them
308   pad1->SetBottomMargin(0);
309   pad2->SetTopMargin(0);
310
311   pad1->cd();
312   gPad->SetGridx();
313   gPad->SetGridy();
314
315   mcHist->GetXaxis()->SetLabelSize(0.06);
316   mcHist->GetYaxis()->SetLabelSize(0.06);
317   mcHist->GetXaxis()->SetTitleSize(0.06);
318   mcHist->GetYaxis()->SetTitleSize(0.06);
319   mcHist->GetYaxis()->SetTitleOffset(0.6);
320
321   mcHist->GetXaxis()->SetRangeUser(0, displayRange);
322
323   mcHist->SetTitle(";true multiplicity;Entries");
324   mcHist->SetStats(kFALSE);
325
326   mcHist->DrawCopy("HIST E");
327   gPad->SetLogy();
328
329   result1->SetLineColor(2);
330   result1->SetMarkerStyle(24);
331   result1->DrawCopy("SAME HISTE");
332
333   result2->SetLineColor(4);
334   result2->SetMarkerColor(4);
335   result2->DrawCopy("SAME HISTE");
336
337   TLegend* legend = new TLegend(0.5, 0.6, 0.95, 0.9);
338   legend->AddEntry(mcHist, "True distribution");
339   legend->AddEntry(result1, "Unfolded distribution (syst)", "P");
340   legend->AddEntry(result2, "Unfolded distribution (normal)", "P");
341   legend->SetFillColor(0);
342   legend->SetTextSize(0.06);
343   legend->Draw();
344
345   pad2->cd();
346   pad2->SetBottomMargin(0.15);
347   //gPad->SetGridx();
348   //gPad->SetGridy();
349
350   result1->GetXaxis()->SetLabelSize(0.06);
351   result1->GetYaxis()->SetLabelSize(0.06);
352   result1->GetXaxis()->SetTitleSize(0.06);
353   result1->GetYaxis()->SetTitleSize(0.06);
354   result1->GetYaxis()->SetTitleOffset(0.6);
355
356   result1->GetXaxis()->SetRangeUser(0, displayRange);
357
358   result1->SetTitle(Form(";%s;Entries", GetMultLabel()));
359   result1->SetStats(kFALSE);
360
361   // calculate ratio
362   result1->Sumw2();
363   TH1* ratio = (TH1*) result1->Clone("ratio");
364   result2->Sumw2();
365   ratio->Divide(ratio, result2, 1, 1, "");
366   ratio->SetLineColor(1);
367   ratio->SetMarkerColor(1);
368   ratio->SetMarkerStyle(0);
369   ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
370   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
371
372   ratio->DrawCopy();
373
374   // get average of ratio
375   Float_t sum = 0;
376   for (Int_t i=2; i<=ratioRange; ++i)
377   {
378     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
379   }
380   sum /= ratioRange-1;
381
382   printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
383
384   TLine* line = new TLine(0, 1, displayRange, 1);
385   line->SetLineWidth(2);
386   line->Draw();
387
388   line = new TLine(0, 1.1, displayRange, 1.1);
389   line->SetLineWidth(2);
390   line->SetLineStyle(2);
391   line->Draw();
392   line = new TLine(0, 0.9, displayRange, 0.9);
393   line->SetLineWidth(2);
394   line->SetLineStyle(2);
395   line->Draw();
396
397   canvas->Modified();
398
399   canvas->SaveAs(epsName);
400
401   return canvas;
402 }
403
404 TCanvas* Draw2ResultRatios(TH1* mcHist, TH1* result1, TH1* result2, TH1* ratio2, TH1* ratio3, TString epsName)
405 {
406   // draws the 3 plots in the upper plot
407   // draws the ratio between result1 and result2 in the lower plot
408   // also draws ratio2 and ratio3 in the lower plot, uses their name for the legend
409
410   // normalize unfolded result to mc hist
411   result1->Scale(1.0 / result1->Integral(2, 200));
412   result1->Scale(mcHist->Integral(2, 200));
413   result2->Scale(1.0 / result2->Integral(2, 200));
414   result2->Scale(mcHist->Integral(2, 200));
415
416   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
417   canvas->Range(0, 0, 1, 1);
418
419   TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
420   pad1->Draw();
421
422   TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
423   pad2->Draw();
424
425   pad1->SetRightMargin(0.05);
426   pad2->SetRightMargin(0.05);
427
428   // no border between them
429   pad1->SetBottomMargin(0);
430   pad2->SetTopMargin(0);
431
432   pad1->cd();
433   gPad->SetGridx();
434   gPad->SetGridy();
435
436   mcHist->GetXaxis()->SetLabelSize(0.06);
437   mcHist->GetYaxis()->SetLabelSize(0.06);
438   mcHist->GetXaxis()->SetTitleSize(0.06);
439   mcHist->GetYaxis()->SetTitleSize(0.06);
440   mcHist->GetYaxis()->SetTitleOffset(0.6);
441
442   mcHist->GetXaxis()->SetRangeUser(0, displayRange);
443
444   mcHist->SetTitle(";True multiplicity;Entries");
445   mcHist->SetStats(kFALSE);
446
447   mcHist->DrawCopy("HIST E");
448   gPad->SetLogy();
449
450   result1->SetLineColor(2);
451   result1->SetMarkerColor(2);
452   result1->SetMarkerStyle(24);
453   result1->DrawCopy("SAME HISTE");
454
455   result2->SetLineColor(4);
456   result2->SetMarkerColor(4);
457   result2->SetMarkerStyle(5);
458   result2->DrawCopy("SAME HISTE");
459
460   TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
461   legend->AddEntry(mcHist, "True distribution");
462   legend->AddEntry(result1, "Unfolded distribution (5%)", "P");
463   legend->AddEntry(result2, "Unfolded distribution (normal)", "P");
464   legend->SetFillColor(0);
465   legend->SetTextSize(0.06);
466   legend->Draw();
467
468   pad2->cd();
469   pad2->SetBottomMargin(0.15);
470   //gPad->SetGridx();
471   //gPad->SetGridy();
472
473   result1->GetXaxis()->SetLabelSize(0.06);
474   result1->GetYaxis()->SetLabelSize(0.06);
475   result1->GetXaxis()->SetTitleSize(0.06);
476   result1->GetYaxis()->SetTitleSize(0.06);
477   result1->GetYaxis()->SetTitleOffset(0.6);
478
479   result1->GetXaxis()->SetRangeUser(0, displayRange);
480
481   result1->SetTitle(Form(";%s;Entries", GetMultLabel()));
482   result1->SetStats(kFALSE);
483
484   // calculate ratio
485   result1->Sumw2();
486   TH1* ratio = (TH1*) result1->Clone("ratio");
487   result2->Sumw2();
488   ratio->Divide(ratio, result2, 1, 1, "");
489   ratio->SetLineColor(1);
490   ratio->SetMarkerColor(1);
491   ratio->SetMarkerStyle(24);
492   ratio->GetYaxis()->SetTitle("Ratio (change / normal)");
493   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
494
495   ratio2->SetLineColor(2);
496   ratio2->SetMarkerColor(2);
497   ratio2->SetMarkerStyle(25);
498
499   ratio3->SetLineColor(4);
500   ratio3->SetMarkerColor(4);
501   ratio3->SetMarkerStyle(26);
502   
503   ratio->DrawCopy();
504   ratio2->DrawCopy("SAME");
505   ratio3->DrawCopy("SAME");
506   
507   legend2 = new TLegend(0.3, 0.8, 0.8, 0.95);
508   legend2->SetNColumns(3);
509   legend2->SetFillColor(0);
510   legend2->SetTextSize(0.06);
511   legend2->AddEntry(ratio, "5% change", "P");
512   legend2->AddEntry(ratio2, "2% change", "P");
513   legend2->AddEntry(ratio3, "1% change", "P");
514   legend2->Draw();
515
516   // get average of ratio
517   Float_t sum = 0;
518   for (Int_t i=2; i<=ratioRange; ++i)
519   {
520     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
521   }
522   sum /= ratioRange-1;
523
524   printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
525
526   TLine* line = new TLine(0, 1, displayRange, 1);
527   line->SetLineWidth(2);
528   line->Draw();
529
530   line = new TLine(0, 1.1, displayRange, 1.1);
531   line->SetLineWidth(2);
532   line->SetLineStyle(2);
533   line->Draw();
534   line = new TLine(0, 0.9, displayRange, 0.9);
535   line->SetLineWidth(2);
536   line->SetLineStyle(2);
537   line->Draw();
538
539   canvas->Modified();
540
541   canvas->SaveAs(epsName);
542
543   return canvas;
544 }
545
546 TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE)
547 {
548   // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
549   // a systematic effect
550
551   // normalize results
552   result->Scale(1.0 / result->Integral(2, 200));
553
554   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 500);
555   canvas->SetTopMargin(0.05);
556   canvas->SetRightMargin(0.05);
557   canvas->SetGridx();
558   canvas->SetGridy();
559
560   result->GetXaxis()->SetRangeUser(0, displayRange);
561   result->GetYaxis()->SetRangeUser(0.55, 1.45);
562   result->SetStats(kFALSE);
563
564   // to get the axis how we want it
565   TH1* dummy = (TH1*) result->Clone("dummy");
566   dummy->Reset();
567   dummy->SetTitle(Form(";%s;Ratio", GetMultLabel()));
568   dummy->DrawCopy();
569   delete dummy;
570
571   Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
572
573   TLegend* legend = new TLegend(0.2, 0.7, 0.7, 0.93);
574   legend->SetFillColor(0);
575   legend->SetTextSize(0.04);
576   if (nResultSyst > 6)
577     legend->SetNColumns(2);
578
579   for (Int_t n=0; n<nResultSyst; ++n)
580   {
581     resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
582
583     // calculate ratio
584     TH1* ratio = (TH1*) result->Clone("ratio");
585     ratio->Divide(ratio, resultSyst[n], 1, 1, "");
586     ratio->GetXaxis()->SetRangeUser(0, displayRange);
587
588     if (firstMarker)
589       ratio->SetMarkerStyle(5);
590
591     ratio->SetLineColor(colors[n / 2]);
592     if ((n % 2))
593       ratio->SetLineStyle(2);
594
595     TString drawStr("SAME HIST");
596     if (n == 0 && firstMarker)
597       drawStr = "SAME P";
598     if (errors)
599       drawStr += " E";
600
601     ratio->DrawCopy(drawStr);
602
603     if (legendStrings && legendStrings[n])
604       legend->AddEntry(ratio, legendStrings[n], "L");
605
606     // get average of ratio
607     Float_t sum = 0;
608     for (Int_t i=2; i<=ratioRange; ++i)
609       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
610     sum /= ratioRange-1;
611
612     printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
613   }
614
615   if (legendStrings)
616     legend->Draw();
617
618   TLine* line = new TLine(-0.5, 1, displayRange, 1);
619   line->SetLineWidth(2);
620   line->Draw();
621
622   line = new TLine(-0.5, 1.1, displayRange, 1.1);
623   line->SetLineWidth(2);
624   line->SetLineStyle(2);
625   line->Draw();
626   line = new TLine(-0.5, 0.9, displayRange, 0.9);
627   line->SetLineWidth(2);
628   line->SetLineStyle(2);
629   line->Draw();
630
631   canvas->SaveAs(epsName);
632   canvas->SaveAs(Form("%s.gif", epsName.Data()));
633
634   return canvas;
635 }
636
637 TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
638 {
639   // draws the ratios of each mc to the corresponding result
640
641   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
642   canvas->SetRightMargin(0.05);
643   canvas->SetTopMargin(0.05);
644
645   for (Int_t n=0; n<nResultSyst; ++n)
646   {
647     // normalize
648     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
649     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
650
651     result[n]->GetXaxis()->SetRangeUser(0, displayRange);
652     result[n]->SetStats(kFALSE);
653
654     // calculate ratio
655     TH1* ratio = (TH1*) result[n]->Clone("ratio");
656     ratio->Divide(mc[n], ratio, 1, 1, "B");
657
658     // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
659     ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
660
661     if (smooth)
662       Smooth(ratio);
663
664     ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
665     ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
666
667     if (dashed)
668     {
669       ratio->SetLineColor((n/2)+1);
670       ratio->SetLineStyle((n%2)+1);
671     }
672     else
673       ratio->SetLineColor(n+1);
674
675     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
676
677     // get average of ratio
678     Float_t sum = 0;
679     for (Int_t i=2; i<=ratioRange; ++i)
680       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
681     sum /= ratioRange-1;
682
683     printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
684   }
685
686   TLine* line = new TLine(0, 1, displayRange, 1);
687   line->SetLineWidth(2);
688   line->Draw();
689
690   line = new TLine(0, 1.1, displayRange, 1.1);
691   line->SetLineWidth(2);
692   line->SetLineStyle(2);
693   line->Draw();
694   line = new TLine(0, 0.9, displayRange, 0.9);
695   line->SetLineWidth(2);
696   line->SetLineStyle(2);
697   line->Draw();
698
699   canvas->Modified();
700
701   canvas->SaveAs(epsName);
702   canvas->SaveAs(Form("%s.gif", epsName.Data()));
703
704   return canvas;
705 }
706
707 TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
708 {
709   // draws the ratios of each mc to the corresponding result
710   // deducts from each ratio the ratio of mcBase / resultBase
711
712   // normalize
713   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
714   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
715
716   // calculate ratio
717   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
718   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
719
720   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
721   canvas->SetRightMargin(0.05);
722   canvas->SetTopMargin(0.05);
723
724   for (Int_t n=0; n<nResultSyst; ++n)
725   {
726     // normalize
727     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
728     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
729
730     result[n]->GetXaxis()->SetRangeUser(0, displayRange);
731     result[n]->SetStats(kFALSE);
732
733     // calculate ratio
734     TH1* ratio = (TH1*) result[n]->Clone("ratio");
735     ratio->Divide(mc[n], ratio, 1, 1, "B");
736     ratio->Add(ratioBase, -1);
737
738     ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
739     ratio->GetYaxis()->SetRangeUser(-1, 1);
740     ratio->SetLineColor(n+1);
741     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
742
743     // get average of ratio
744     Float_t sum = 0;
745     for (Int_t i=2; i<=ratioRange; ++i)
746       sum += TMath::Abs(ratio->GetBinContent(i));
747     sum /= ratioRange-1;
748
749     printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, ratioRange, sum);
750   }
751
752   TLine* line = new TLine(0, 0, displayRange, 0);
753   line->SetLineWidth(2);
754   line->Draw();
755
756   line = new TLine(0, 0.1, displayRange, 0.1);
757   line->SetLineWidth(2);
758   line->SetLineStyle(2);
759   line->Draw();
760   line = new TLine(0, -0.1, displayRange, -0.1);
761   line->SetLineWidth(2);
762   line->SetLineStyle(2);
763   line->Draw();
764
765   canvas->Modified();
766
767   canvas->SaveAs(epsName);
768   canvas->SaveAs(Form("%s.gif", epsName.Data()));
769
770   return canvas;
771 }
772
773 TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
774 {
775   // draws the ratios of each mc to the corresponding result
776   // deducts from each ratio the ratio of mcBase / resultBase
777   // smoothens the ratios by a sliding window
778
779   // normalize
780   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
781   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
782
783   // calculate ratio
784   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
785   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
786
787   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
788   canvas->SetRightMargin(0.05);
789   canvas->SetTopMargin(0.05);
790
791   for (Int_t n=0; n<nResultSyst; ++n)
792   {
793     // normalize
794     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
795     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
796
797     result[n]->GetXaxis()->SetRangeUser(0, displayRange);
798     result[n]->SetStats(kFALSE);
799
800     // calculate ratio
801     TH1* ratio = (TH1*) result[n]->Clone("ratio");
802     ratio->Divide(mc[n], ratio, 1, 1, "B");
803     ratio->Add(ratioBase, -1);
804
805     //new TCanvas; ratio->DrawCopy();
806     // clear 0 bin
807     ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
808
809     Smooth(ratio);
810
811     //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
812
813     canvas->cd();
814     ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
815     ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
816     ratio->SetLineColor((n / 2)+1);
817     ratio->SetLineStyle((n % 2)+1);
818     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
819
820     // get average of ratio
821     Float_t sum = 0;
822     for (Int_t i=2; i<=150; ++i)
823       sum += TMath::Abs(ratio->GetBinContent(i));
824     sum /= 149;
825
826     printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
827   }
828
829   TLine* line = new TLine(0, 0, displayRange, 0);
830   line->SetLineWidth(2);
831   line->Draw();
832
833   line = new TLine(0, 0.1, displayRange, 0.1);
834   line->SetLineWidth(2);
835   line->SetLineStyle(2);
836   line->Draw();
837   line = new TLine(0, -0.1, displayRange, -0.1);
838   line->SetLineWidth(2);
839   line->SetLineStyle(2);
840   line->Draw();
841
842   canvas->Modified();
843
844   canvas->SaveAs(epsName);
845   canvas->SaveAs(Form("%s.gif", epsName.Data()));
846
847   return canvas;
848 }
849
850 void DrawResiduals(const char* fileName, const char* epsName)
851 {
852   loadlibs();
853
854   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
855   TFile::Open(fileName);
856   mult->LoadHistograms("Multiplicity");
857
858   TH1* measured = mult->GetMultiplicityESD(etaRange)->ProjectionY("myesd", 1, 1);
859   TH1* unfoldedFolded = mult->CalculateMultiplicityESD(mult->GetMultiplicityESDCorrected(etaRange), etaRange)->ProjectionY("myfolded", 1, 1);
860   
861   // normalize
862   unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, displayRange+1));
863   unfoldedFolded->Scale(measured->Integral(2, displayRange+1));
864
865   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
866   canvas->Range(0, 0, 1, 1);
867
868   TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 1, 1);
869   pad1->Draw();
870   pad1->SetGridx();
871   pad1->SetGridy();
872
873   TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 1, 0.5);
874   pad2->Draw();
875   pad2->SetGridx();
876   pad2->SetGridy();
877
878   TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
879   pad3->SetGridx();
880   pad3->SetGridy();
881   pad3->SetRightMargin(0.05);
882   pad3->SetTopMargin(0.05);
883   pad3->Draw();
884
885   pad1->SetRightMargin(0.05);
886   pad2->SetRightMargin(0.05);
887
888   // no border between them
889   pad1->SetBottomMargin(0);
890   pad2->SetTopMargin(0);
891
892   pad1->cd();
893
894   measured->GetXaxis()->SetLabelSize(0.06);
895   measured->GetYaxis()->SetLabelSize(0.06);
896   measured->GetXaxis()->SetTitleSize(0.06);
897   measured->GetYaxis()->SetTitleSize(0.06);
898   measured->GetYaxis()->SetTitleOffset(0.6);
899
900   measured->GetXaxis()->SetRangeUser(0, displayRange);
901
902   measured->SetTitle(Form(";%s;Entries", GetMultLabel(etaRange, kFALSE)));
903   measured->SetStats(kFALSE);
904
905   measured->DrawCopy("HIST");
906   gPad->SetLogy();
907
908   unfoldedFolded->SetMarkerStyle(5);
909   unfoldedFolded->SetMarkerColor(2);
910   unfoldedFolded->SetLineColor(2);
911   unfoldedFolded->DrawCopy("SAME PHIST");
912
913   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
914   legend->AddEntry(measured, "Measured distribution", "L");
915   legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution", "P");
916   legend->SetFillColor(0);
917   legend->SetTextSize(0.06);
918   legend->Draw();
919
920   pad2->cd();
921   pad2->SetBottomMargin(0.15);
922
923   // calculate ratio
924   measured->Sumw2();
925   TH1* residual = (TH1*) measured->Clone("residual");
926   unfoldedFolded->Sumw2();
927
928   residual->Add(unfoldedFolded, -1);
929
930   // projection
931   TH1* residualHist = new TH1F("residualHist", ";", 11, -3, 3);
932   residualHist->Sumw2();
933
934   Float_t chi2 = 0;
935   for (Int_t i=1; i<=displayRange+1; ++i)
936   {
937     if (measured->GetBinError(i) > 0)
938     {
939       residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
940       residual->SetBinError(i, 1);
941
942       residualHist->Fill(residual->GetBinContent(i));
943       chi2 += residual->GetBinContent(i) * residual->GetBinContent(i);
944     }
945     else
946     {
947       residual->SetBinContent(i, 0);
948       residual->SetBinError(i, 0);
949     }
950   }
951   
952   Printf("chi2 / ndf = %f / %d = %f", chi2, displayRange+1, chi2 / (displayRange+1));
953
954   residual->GetYaxis()->SetTitle("Residuals:   (1/e) (M - R  #otimes U)");
955   residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
956   residual->DrawCopy();
957
958   TLine* line = new TLine(-0.5, 0, displayRange + 0.5, 0);
959   line->SetLineWidth(2);
960   line->Draw();
961
962   pad3->cd();
963   residualHist->SetStats(kFALSE);
964   residualHist->SetLabelSize(0.08, "xy");
965   residualHist->Fit("gaus");
966   residualHist->Draw("HIST");
967   residualHist->FindObject("gaus")->Draw("SAME");
968
969   canvas->Modified();
970   canvas->SaveAs(canvas->GetName());
971
972   //const char* epsName2 = "proj.eps";
973   //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
974   //canvas->SetGridx();
975   //canvas->SetGridy();
976
977   //canvas->SaveAs(canvas->GetName());
978 }
979
980 void chi2FluctuationResult()
981 {
982   loadlibs();
983
984   TFile::Open("chi2_noregularization.root");
985   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
986   mult->LoadHistograms("Multiplicity");
987
988   TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
989
990   mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
991
992   TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_1");
993   ((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetRangeUser(0, displayRange);
994   ((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->SetRangeUser(0, displayRange);
995   //((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetTitle(GetMultTitle());
996   //((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->->SetTitle(GetMultTitle(etaRange, kFALSE));
997   canvas->SaveAs("chi2FluctuationResult.eps");
998 }
999
1000 void DrawUnfolded(const char* fileName, const char* eps)
1001 {
1002   loadlibs();
1003   
1004   TFile::Open(fileName);
1005
1006   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1007   mult->LoadHistograms("Multiplicity");
1008
1009   TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult->GetMultiplicityVtx(etaRange)->GetNbinsX());
1010   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1011
1012   DrawResultRatio(mcHist, result, eps);
1013 }
1014
1015 void minimizationInfluenceAlpha()
1016 {
1017   loadlibs();
1018
1019   TFile::Open(measuredFile);
1020   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1021   mult2->LoadHistograms("Multiplicity");
1022
1023   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult2->GetMultiplicityVtx(etaRange)->GetNbinsX());
1024   mcHist->Scale(1.0 / mcHist->Integral());
1025   mcHist->SetStats(kFALSE);
1026   mcHist->SetTitle(";True multiplicity n in |#eta| < 1.0;P(n)");
1027
1028   TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 350);
1029   canvas->Divide(3, 1, 0.005);
1030
1031   TFile::Open("chi2compare-influencealpha/EvaluateChi2Method.root");
1032   
1033   TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_3.162278");
1034   TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_07_2_10000.000000");
1035   TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_13_2_10000000.000000");
1036
1037   /*mcHist->Rebin(2);  mcHist->Scale(0.5);
1038   hist1->Rebin(2);   hist1->Scale(0.5);
1039   hist2->Rebin(2);   hist2->Scale(0.5);
1040   hist3->Rebin(2);   hist3->Scale(0.5);*/
1041
1042   mcHist->GetXaxis()->SetRangeUser(0, displayRange);
1043   mcHist->SetLabelSize(0.06, "xy");
1044   mcHist->SetTitleSize(0.06, "xy");
1045   mcHist->GetYaxis()->SetTitleOffset(1.5);
1046   
1047   canvas->cd(1);
1048   
1049   gPad->SetLogy();
1050   gPad->SetRightMargin(0.03);
1051   gPad->SetLeftMargin(0.19);
1052   gPad->SetTopMargin(0.05);
1053   gPad->SetBottomMargin(0.13);
1054   gPad->SetGridx();
1055   gPad->SetGridy();
1056   mcHist->Draw();
1057   hist1->SetMarkerStyle(5);
1058   hist1->SetMarkerColor(2);
1059   hist1->Draw("SAME PE");
1060
1061   canvas->cd(2);
1062   gPad->SetRightMargin(0.03);
1063   gPad->SetLeftMargin(0.19);
1064   gPad->SetTopMargin(0.05);
1065   gPad->SetBottomMargin(0.13);
1066   gPad->SetLogy();
1067   gPad->SetGridx();
1068   gPad->SetGridy();
1069   mcHist->Draw();
1070   hist2->SetMarkerStyle(5);
1071   hist2->SetMarkerColor(2);
1072   hist2->Draw("SAME PE");
1073
1074   canvas->cd(3);
1075   gPad->SetRightMargin(0.03);
1076   gPad->SetLeftMargin(0.19);
1077   gPad->SetTopMargin(0.05);
1078   gPad->SetBottomMargin(0.13);
1079   gPad->SetLogy();
1080   gPad->SetGridx();
1081   gPad->SetGridy();
1082   mcHist->Draw();
1083   hist3->SetMarkerStyle(5);
1084   hist3->SetMarkerColor(2);
1085   hist3->Draw("SAME PE");
1086
1087   canvas->SaveAs("minimizationInfluenceAlpha.eps");
1088 }
1089
1090 void NBDFit()
1091 {
1092   gSystem->Load("libPWG0base");
1093
1094   TFile::Open(correctionFile);
1095   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1096   mult->LoadHistograms("Multiplicity");
1097
1098   TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
1099   fCurrentESD->Sumw2();
1100   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1101
1102   TF1* 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])");
1103   func->SetParNames("scaling", "averagen", "k");
1104   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
1105   func->SetParLimits(1, 0.001, 1000);
1106   func->SetParLimits(2, 0.001, 1000);
1107   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
1108
1109   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
1110   lognormal->SetParNames("scaling", "mean", "sigma");
1111   lognormal->SetParameters(1, 1, 1);
1112   lognormal->SetParLimits(0, 0, 10);
1113   lognormal->SetParLimits(1, 0, 100);
1114   lognormal->SetParLimits(2, 1e-3, 10);
1115
1116   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
1117   fCurrentESD->SetStats(kFALSE);
1118   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
1119   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
1120   fCurrentESD->Draw("HIST");
1121   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1122   fCurrentESD->Fit(func, "W0", "", 0, 50);
1123   func->SetRange(0, 100);
1124   func->Draw("SAME");
1125   printf("chi2 = %f\n", func->GetChisquare());
1126
1127   fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
1128   lognormal->SetLineColor(2);
1129   lognormal->SetLineStyle(2);
1130   lognormal->SetRange(0, 100);
1131   lognormal->Draw("SAME");
1132
1133   canvas->SaveAs("NBDFit.eps");
1134 }
1135
1136 void StartingConditions()
1137 {
1138   // data generated by runMultiplicitySelector.C StartingConditions
1139
1140   const char* name = "StartingConditions";
1141
1142   TFile* file = TFile::Open(Form("%s.root", name));
1143
1144   TCanvas* canvas = new TCanvas(name, name, 1200, 600);
1145   canvas->Divide(2, 1);
1146
1147   TH1* mc = (TH1*) file->Get("mc");
1148   mc->Sumw2();
1149   mc->Scale(1.0 / mc->Integral(2, displayRange));
1150
1151   //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
1152
1153   TLegend* legend = new TLegend(0.6, 0.7, 0.99, 0.99);
1154   legend->SetFillColor(0);
1155   legend->SetTextSize(0.04);
1156   
1157   const char* names[] = { "True", "Measured 1", "Measured 2", "Measured 3", "NBD", "Flat" };
1158
1159   for (Int_t i=0; i<6; ++i)
1160   {
1161     Int_t id = i;
1162     if (id > 2)
1163       id += 2;
1164
1165     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1166     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
1167     
1168     chi2Result->Scale(1.0 / chi2Result->Integral(2, displayRange));
1169     bayesResult->Scale(1.0 / bayesResult->Integral(2, displayRange));
1170
1171     chi2Result->Divide(mc, chi2Result, 1, 1, "");
1172     bayesResult->Divide(mc, bayesResult, 1, 1, "");
1173
1174     chi2Result->SetTitle(Form("a) #chi^{2}-minimization;%s;MC / unfolded", GetMultLabel()));
1175     chi2Result->GetXaxis()->SetRangeUser(0, displayRange);
1176     chi2Result->GetYaxis()->SetRangeUser(0.7, 1.3);
1177     chi2Result->GetYaxis()->SetTitleOffset(1.7);
1178     //chi2Result->SetMarkerStyle(marker[i]);
1179     chi2Result->SetLineColor(i+1);
1180     chi2Result->SetMarkerColor(i+1);
1181     chi2Result->SetStats(kFALSE);
1182
1183     bayesResult->SetTitle(Form("b) Bayesian unfolding;%s;MC / unfolded", GetMultLabel()));
1184     bayesResult->GetXaxis()->SetRangeUser(0, displayRange);
1185     bayesResult->GetYaxis()->SetRangeUser(0.7, 1.3);
1186     bayesResult->GetYaxis()->SetTitleOffset(1.7);
1187     bayesResult->SetStats(kFALSE);
1188     //bayesResult->SetLineColor(2);
1189     bayesResult->SetLineColor(i+1);
1190
1191     canvas->cd(1);
1192     gPad->SetRightMargin(0.05);
1193     gPad->SetLeftMargin(0.12);
1194     gPad->SetGridx();
1195     gPad->SetGridy();
1196     chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1197
1198     canvas->cd(2);
1199     gPad->SetRightMargin(0.05);
1200     gPad->SetLeftMargin(0.12);
1201     gPad->SetGridx();
1202     gPad->SetGridy();
1203     bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1204
1205     //TLine* line = new TLine(0, 1, 150, 1);
1206     //line->Draw();
1207
1208     legend->AddEntry(chi2Result, names[i]);
1209   }
1210
1211   canvas->cd(1);
1212   legend->Draw();
1213   canvas->cd(2);
1214   legend->Draw();
1215   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1216 }
1217
1218 void StatisticsPlot()
1219 {
1220   const char* name = "StatisticsPlot";
1221
1222   TFile* file = TFile::Open(Form("%s.root", name));
1223
1224   TCanvas* canvas = new TCanvas(name, name, 600, 400);
1225
1226   TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
1227   fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1228   fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1229   fitResultsChi2->Draw("AP");
1230
1231   TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1232   fitResultsChi2->Fit(f);
1233
1234   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1235
1236   TH1* mc[5];
1237   TH1* result[5];
1238
1239   const char* plotname = "chi2Result";
1240
1241   name = "StatisticsPlotRatios";
1242   canvas = new TCanvas(name, name, 600, 400);
1243
1244   for (Int_t i=0; i<5; ++i)
1245   {
1246     mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1247     result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1248
1249     result[i]->SetLineColor(i+1);
1250     result[i]->Draw(((i == 0) ? "" : "SAME"));
1251   }
1252 }
1253
1254 void Draw2Unfolded(const char* file1, const char* file2, const char* output)
1255 {
1256   loadlibs();
1257   
1258   TFile::Open(file1);
1259   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1260   mult->LoadHistograms("Multiplicity");
1261
1262   // result with systematic effect
1263   TFile::Open(file2);
1264   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1265   mult2->LoadHistograms("Multiplicity");
1266
1267   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1268   TH1* result1 = (TH1*) mult2->GetMultiplicityESDCorrected(etaRange)->Clone("result1"); // from file2 (with syst)
1269   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");  // from file1 (without syst)
1270
1271   DrawResultRatio(mcHist, result1, "tmp1.eps");
1272   DrawResultRatio(mcHist, result2, "tmp2.eps");
1273   Draw2ResultRatio(mcHist, result1, result2, output);
1274 }
1275
1276 void PythiaPhojet()
1277 {
1278   loadlibs();
1279   
1280   displayRange = 55;
1281   Draw2Unfolded("self.root", "pythia.root", "test.eps");
1282   
1283   canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last();
1284   pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First();
1285   pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last();
1286   legend = (TLegend*)pad1->GetListOfPrimitives()->Last();
1287   
1288   ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (Pythia / Phojet)");
1289   ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (Pythia)");
1290   ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (Phojet)");
1291   canvas->SaveAs("PythiaPhojet.eps");
1292 }
1293
1294 void Misalignment()
1295 {
1296   loadlibs();
1297   
1298   Draw2Unfolded("chi2_ideal.root", "chi2_misaligned.root", "test.eps");
1299   
1300   canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last();
1301   pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First();
1302   pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last();
1303   legend = (TLegend*)pad1->GetListOfPrimitives()->Last();
1304
1305   ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (misaligned / realigned)");
1306   ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (misaligned)");
1307   ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (realigned)");
1308   canvas->SaveAs("SystematicMisalignment.eps");
1309 }
1310
1311 void SystematicLowEfficiency()
1312 {
1313   Draw2Unfolded("chi2.root", "chi2_loweff_5.root", "SystematicLowEfficiency.eps");
1314 }
1315
1316 void SystematicLowEfficiency2()
1317 {
1318   loadlibs();
1319   
1320   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("chi2.root");
1321   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1322   TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1323   result2->Scale(1.0 / result2->Integral(2, displayRange));
1324   result2->Scale(mcHist->Integral(2, displayRange));
1325   
1326   mult = AliMultiplicityCorrection::Open("chi2_loweff_5.root");
1327   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1328   
1329   mult = AliMultiplicityCorrection::Open("chi2_loweff_2.root");
1330   TH1* ratio2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio2");
1331   ratio2->Scale(1.0 / ratio2->Integral(2, displayRange));
1332   ratio2->Scale(mcHist->Integral(2, displayRange));
1333   ratio2->Divide(result2);
1334   
1335   mult = AliMultiplicityCorrection::Open("chi2_loweff_1.root");
1336   TH1* ratio3 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio3");
1337   ratio3->Scale(1.0 / ratio3->Integral(2, displayRange));
1338   ratio3->Scale(mcHist->Integral(2, displayRange));
1339   ratio3->Divide(result2);
1340   
1341   Draw2ResultRatios(mcHist, result1, result2, ratio2, ratio3, "SystematicLowEfficiency2.eps");
1342 }
1343
1344 void SystematicMisalignment()
1345 {
1346   gSystem->Load("libPWG0base");
1347
1348   TFile::Open(correctionFile);
1349   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1350   mult->LoadHistograms("Multiplicity");
1351
1352   TFile::Open("multiplicityMC_100k_fullmis.root");
1353   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1354   mult2->LoadHistograms("Multiplicity");
1355
1356   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1357   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1358   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1359
1360   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1361   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1362
1363   DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1364 }
1365
1366 void SystematicMisalignmentTPC()
1367 {
1368   gSystem->Load("libPWG0base");
1369
1370   SetTPC();
1371
1372   TFile::Open(correctionFile);
1373   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1374   mult->LoadHistograms("Multiplicity");
1375
1376   TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1377   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1378   mult2->LoadHistograms("Multiplicity");
1379
1380   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1381   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1382   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1383
1384   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1385   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1386
1387   DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1388 }
1389
1390 void LowMomentumEffectSPD()
1391 {
1392   // this function increases/reduces the correction as function of pt between 0 and 0.2 gev by +-50% to 0% and checks the effect on the overall correction factor
1393   // only a normal acceptance region is considered to not get a bias by the edges
1394   
1395   loadlibs();
1396   TFile::Open("multiplicity.root");
1397   
1398     
1399   AliCorrection* correction[8];
1400   Float_t values[3];
1401   
1402   for (Int_t loop=0; loop<3; loop++)
1403   {
1404     Float_t sumGen = 0;
1405     Float_t sumMeas = 0;
1406     
1407     Printf("loop %d", loop);
1408     for (Int_t i=0; i<4; ++i)
1409     {
1410       Printf("correction %d", i);
1411   
1412       TString name; name.Form("correction_%d", i);
1413       correction[i] = new AliCorrection(name, name);
1414       correction[i]->LoadHistograms();
1415       
1416       TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1417       TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1418   
1419       Float_t vtxRange = 5.9;
1420       gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1421       meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1422       
1423       Float_t etaRange = 0.99;
1424       gene->GetYaxis()->SetRangeUser(-etaRange, etaRange);
1425       meas->GetYaxis()->SetRangeUser(-etaRange, etaRange);
1426   
1427       TH1* genePt = gene->Project3D(Form("z_%d", i));
1428       TH1* measPt = meas->Project3D(Form("z_%d", i));
1429   
1430       if (loop > 0)
1431       {
1432         for (Int_t x=1; x<=genePt->GetNbinsX(); x++)
1433         {
1434           Float_t pt = genePt->GetXaxis()->GetBinCenter(x);
1435           //Printf("%f", pt);
1436           if (pt < 0.2)
1437           {
1438             Float_t factor = 1;
1439             if (loop == 1)
1440               factor = 1.5 - pt / 0.2 * 0.5;
1441             if (loop == 2)
1442               factor = 0.5 + pt / 0.2 * 0.5;
1443             //Printf("%f", factor);
1444             genePt->SetBinContent(x, genePt->GetBinContent(x) * factor);
1445             measPt->SetBinContent(x, measPt->GetBinContent(x) * factor);
1446           }
1447         }
1448       }
1449       
1450       //new TCanvas; genePt->DrawCopy(); measPt->DrawCopy("SAME");
1451   
1452       sumGen += genePt->Integral();
1453       sumMeas += measPt->Integral();  
1454       
1455       Float_t average = measPt->Integral() / genePt->Integral();
1456       
1457       Printf("The average efficiency of this correction is %f", average);
1458     }
1459     
1460     Float_t average = sumMeas / sumGen;
1461       
1462     Printf("The average efficiency of all corrections is %f", average);
1463     values[loop] = average;
1464   }
1465   
1466   Printf("relative is %f and %f", values[1] / values[0], values[2] / values[0]);
1467 }
1468   
1469
1470 void EfficiencySpecies(Bool_t addDecayStopped = kFALSE)
1471 {
1472   loadlibs();
1473
1474   Int_t marker[] = {24, 25, 26, 27};
1475   Int_t color[] = {1, 2, 4, 3};
1476
1477   // SPD TPC
1478   //const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
1479   //const char* fileName[] = { "spd/multiplicity.root", "tpc/multiplicity.root" };
1480   const char* fileName[] = { "multiplicity.root", "multiplicity.root" };
1481   Float_t etaRangeArr[] = {0.49, 0.9};
1482   const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1483
1484   TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1485   canvas->Divide(2, 1);
1486
1487   TCanvas* canvas3 = new TCanvas("EfficiencySpecies_comb", "EfficiencySpecies_comb", 600, 600);
1488   gPad->SetGridx();
1489   gPad->SetGridy();
1490   gPad->SetRightMargin(0.05);
1491   gPad->SetTopMargin(0.05);
1492   
1493   TLegend* legends[2];
1494   
1495   for (Int_t loop=1; loop<2; ++loop)
1496   {
1497     Printf("%s", fileName[loop]);
1498
1499     TCanvas* canvas2 = new TCanvas(Form("EfficiencySpecies_%d", loop), Form("EfficiencySpecies_%d", loop), 600, 600);
1500     gPad->SetGridx();
1501     gPad->SetGridy();
1502     gPad->SetRightMargin(0.05);
1503     gPad->SetTopMargin(0.05);
1504     
1505     AliCorrection* correction[8];
1506
1507     canvas->cd(loop+1);
1508
1509     gPad->SetGridx();
1510     gPad->SetGridy();
1511     gPad->SetRightMargin(0.05);
1512
1513     TLegend* legend = new TLegend(0.6, 0.4, 0.85, 0.6);
1514     legend->SetFillColor(0);
1515     legend->SetEntrySeparation(0.2);
1516     legend->SetTextSize(gStyle->GetTextSize());
1517     
1518     legends[loop] = new TLegend(0.4+loop*0.3, 0.2, 0.6+loop*0.3, 0.5);
1519     legends[loop]->SetFillColor(0);
1520     legends[loop]->SetEntrySeparation(0.2);
1521     legends[loop]->SetTextSize(gStyle->GetTextSize());
1522     legends[loop]->SetHeader((loop == 0) ? "SPD" : "TPC");
1523
1524     Float_t below = 0;
1525     Float_t total = 0;
1526
1527     TFile* file = TFile::Open(fileName[loop]);
1528     if (!file)
1529     {
1530       Printf("Could not open %s", fileName[loop]);
1531       return;
1532     }
1533
1534     Float_t sumGen = 0;
1535     Float_t sumMeas = 0;
1536
1537     for (Int_t i=0; i<3; ++i)
1538     {
1539       Printf("correction %d", i);
1540
1541       TString name; name.Form("correction_%d", i);
1542       correction[i] = new AliCorrection(name, name);
1543       correction[i]->LoadHistograms();
1544       
1545       TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1546       TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1547
1548       // limit vtx axis
1549       Float_t vtxRange = 3.9;
1550       gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1551       meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1552
1553       // empty over/underflow bin in eta, setting range to +-2 is not enough because this is the maximum range, Project3D takes them into account then (might be a bug)
1554       /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1555         for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1556         {
1557           gene->SetBinContent(x, 0, z, 0);
1558           gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1559           meas->SetBinContent(x, 0, z, 0);
1560           meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1561         }*/
1562
1563       // limit eta axis
1564       Float_t etaBegin = -etaRangeArr[loop];
1565       Float_t etaEnd   = etaRangeArr[loop];
1566       //etaBegin = 0.01;
1567       //etaEnd = -0.01;
1568       gene->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
1569       meas->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
1570
1571       TH1* genePt = gene->Project3D(Form("z_%d", i));
1572       TH1* measPt = meas->Project3D(Form("z_%d", i));
1573
1574       genePt->Sumw2();
1575       measPt->Sumw2();
1576       
1577       for (Int_t x=0; x<=genePt->GetNbinsX()+1; x++)
1578       {
1579         genePt->SetBinError(x, TMath::Sqrt(genePt->GetBinContent(x)));
1580         measPt->SetBinError(x, TMath::Sqrt(measPt->GetBinContent(x)));
1581       }
1582       
1583       sumGen += genePt->Integral();
1584       sumMeas += measPt->Integral();
1585
1586       TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1587       effPt->Reset();
1588       effPt->Divide(measPt, genePt, 1, 1, "B");
1589
1590       Int_t bin = 0;
1591       for (bin=20; bin>=1; bin--)
1592       {
1593         if (effPt->GetBinContent(bin) < 0.5)
1594           break;
1595       }
1596
1597       Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1598
1599       Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1600       Printf("%.4f of the particles are below that momentum", fraction);
1601
1602       below += genePt->Integral(1, bin);
1603       total += genePt->Integral();
1604       
1605       effPt->SetLineColor(color[i]);
1606       effPt->SetMarkerColor(color[i]);
1607       effPt->SetMarkerStyle(marker[i]);
1608       effPt->SetMarkerSize(2);
1609
1610       effPt->GetXaxis()->SetRangeUser(0, 1);
1611       effPt->GetYaxis()->SetRangeUser(0.001, 1);
1612
1613       effPt->GetXaxis()->SetTitleOffset(1.1);
1614       effPt->GetYaxis()->SetTitleOffset(1.2);
1615       
1616       effPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1617
1618       effPt->SetStats(kFALSE);
1619       effPt->SetTitle(titles[loop]);
1620       effPt->GetYaxis()->SetTitle("Efficiency");
1621
1622       canvas->cd(loop+1);
1623       effPt->DrawCopy((i == 0) ? "" : "SAME");
1624  
1625       canvas2->cd();
1626       effPt->SetTitle("");
1627       effPt->DrawCopy((i == 0) ? "" : "SAME");
1628       
1629       canvas3->cd();
1630       effPtClone = (TH1*) effPt->Clone("effPtClone");
1631       effPtClone->SetMarkerStyle(marker[i]-4*loop);
1632       effPtClone->DrawCopy((i == 0 && loop == 0) ? "" : "SAME");
1633
1634       legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1635       legends[loop]->AddEntry(effPtClone, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1636       //legend2->AddEntry(effPt, Form("%s %s", (loop == 0) ? "SPD" : "TPC", ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}"))), "P");
1637
1638       if (addDecayStopped)
1639       {
1640         name.Form("correction_%d", i+4);
1641         corr = new AliCorrection(name, name);
1642         corr->LoadHistograms();
1643         
1644         TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
1645         TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
1646         
1647         // limit axes
1648         gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1649         meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1650         gene->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1651         meas->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1652         
1653         TH1* decayed = gene->Project3D(Form("z_%d", i+4));
1654         TH1* stopped = meas->Project3D(Form("z_%d", i+4));
1655         
1656         Printf("%d: %d decayed, %d stopped, out of %d", i, (Int_t) decayed->Integral(), (Int_t) stopped->Integral(), (Int_t) genePt->Integral());
1657         
1658         decayed->Divide(decayed, genePt, 1, 1, "B");
1659         stopped->Divide(stopped, genePt, 1, 1, "B");
1660         
1661         decayed->SetMarkerStyle(20);
1662         stopped->SetMarkerStyle(21);
1663         stopped->SetMarkerColor(2);
1664         
1665         new TCanvas(Form("all_%d_%d", loop, i), Form("all_%d_%d", loop, i), 600, 600);
1666         effPt->DrawCopy();
1667         decayed->DrawCopy("SAME");
1668         stopped->DrawCopy("SAME");
1669         
1670         decayed->Add(stopped);
1671         decayed->Add(effPt);
1672         decayed->SetMarkerStyle(22);
1673         decayed->SetMarkerColor(4);
1674         decayed->DrawCopy("SAME");
1675       }
1676       
1677     }
1678
1679     Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1680
1681     Printf("%f measured, %f generated, effiency: %f", sumGen, sumMeas, sumMeas / sumGen);
1682
1683     canvas->cd(loop+1);
1684     legend->Draw();
1685   
1686     canvas2->cd();
1687     legend->Draw();
1688     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));  
1689   }
1690
1691   canvas->cd();
1692   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1693
1694   canvas3->cd();
1695   legends[0]->Draw();
1696   legends[1]->Draw();
1697   canvas3->SaveAs(Form("%s.eps", canvas3->GetName()));
1698 }
1699
1700 void DrawpTCutOff()
1701 {
1702 /*
1703 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt_ref.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1704 mv unfolded.root chi2_ptref.root
1705 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1706 mv unfolded.root chi2_pt0.root
1707 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1708 mv unfolded.root chi2_pt1.root
1709 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1710 mv unfolded.root chi2_pt0_25.root
1711 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1712 mv unfolded.root chi2_pt1_25.root
1713 */
1714
1715   loadlibs();
1716   
1717   TH1* results[10];
1718   const char* files[] = { "chi2_ptref.root", "chi2_pt0.root", "chi2_pt1.root", "chi2_pt0_25.root", "chi2_pt1_25.root"};
1719   
1720   Int_t nMax = 5;
1721   for (Int_t i = 0; i<nMax; ++i)
1722   {
1723     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(files[i]);
1724     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1725   }
1726   
1727   const char* legendStrings[] = { "Reduced 50%", "Enhanced 50%", "Reduced 25%", "Enhanced 25%" };
1728   DrawRatio(results[0], nMax-1, results+1, "LowMomentumSyst.eps", kFALSE, legendStrings);
1729 }
1730
1731 void ParticleSpeciesComparison()
1732 {
1733   loadlibs();
1734
1735   TH1* results[10];
1736   TH1* mc = 0;
1737   
1738   // loop over cases (normal, enhanced/reduced ratios)
1739   Int_t nMax = 9;
1740   for (Int_t i = 0; i<nMax; ++i)
1741   {
1742     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(Form("chi2_species_%d.root", i), Form("Multiplicity_%d", i));
1743     if (i == 0)
1744       mc = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymchist", 1, 1);
1745     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1746   }
1747
1748   DrawResultRatio(mc, results[0], "ParticleSpeciesComparison1_1.eps");
1749
1750   for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1751   {
1752     results[0]->SetBinError(i, 0);
1753     mc->SetBinError(i, 0);
1754   }
1755
1756   const char* legendStrings[] = { "K #times 0.5", "K #times 1.5", "p #times 0.5", "p #times 1.5", "K #times 0.5, p #times 0.5", "K #times 1.5, p #times 1.5", "K #times 0.5, p #times 1.5", "K #times 1.5, p #times 0.5" };
1757
1758   DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1759
1760   //not valid: draw chi2 uncertainty on top!
1761   /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1762   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1763   errorHist->SetLineColor(1);
1764   errorHist->SetLineWidth(2);
1765   TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1766   for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1767   {
1768     errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1769     errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1770   }
1771
1772   errorHist->DrawCopy("SAME");
1773   errorHist2->DrawCopy("SAME");*/
1774
1775   //canvas->SaveAs(canvas->GetName());
1776
1777   DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1778
1779   //errorHist->DrawCopy("SAME");
1780   //errorHist2->DrawCopy("SAME");
1781
1782   //canvas2->SaveAs(canvas2->GetName());
1783 }
1784
1785 /*void ParticleSpeciesComparison2()
1786 {
1787   gSystem->Load("libPWG0base");
1788
1789   const char* fileNameMC = "multiplicityMC_400k_syst.root";
1790   const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1791   Bool_t chi2 = 0;
1792
1793   TFile::Open(fileNameMC);
1794   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1795   mult->LoadHistograms();
1796
1797   TH1* mc[10];
1798   TH1* results[10];
1799
1800   // loop over cases (normal, enhanced/reduced ratios)
1801   Int_t nMax = 7;
1802   for (Int_t i = 0; i<nMax; ++i)
1803   {
1804     TString folder;
1805     folder.Form("Multiplicity_%d", i);
1806
1807     AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1808
1809     TFile::Open(fileNameESD);
1810     mult2->LoadHistograms();
1811
1812     mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1813
1814     if (chi2)
1815     {
1816       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1817       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1818       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1819     }
1820     else
1821     {
1822       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1823       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1824     }
1825
1826     //Float_t averageRatio = 0;
1827     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1828
1829     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1830
1831     TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1832     mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1833
1834     //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1835     //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1836
1837     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1838   }
1839
1840   DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1841 }*/
1842
1843 TH1* Invert(TH1* eff)
1844 {
1845   // calculate corr = 1 / eff
1846
1847   TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1848   corr->Reset();
1849
1850   for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1851   {
1852     if (eff->GetBinContent(i) > 0)
1853     {
1854       corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1855       corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1856     }
1857   }
1858
1859   return corr;
1860 }
1861
1862 void TriggerVertexCorrection()
1863 {
1864   //
1865   // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1866   //
1867
1868   loadlibs();
1869
1870   TFile::Open(correctionFile);
1871   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1872   mult->LoadHistograms("Multiplicity");
1873
1874   TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1875   TH1* corrNSD = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kNSD));
1876   TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1877
1878   TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 500);
1879   gPad->SetGridx();
1880   gPad->SetGridy();
1881   gPad->SetTopMargin(0.05);
1882   gPad->SetRightMargin(0.05);
1883
1884   corrINEL->SetStats(kFALSE);
1885   corrINEL->GetXaxis()->SetRangeUser(0, 12);
1886   corrINEL->GetYaxis()->SetRangeUser(0, 8);
1887   corrINEL->SetTitle(Form(";%s;Correction factor", GetMultLabel()));
1888   corrINEL->SetMarkerStyle(22);
1889   corrINEL->Draw("PE");
1890
1891   corrMB->SetStats(kFALSE);
1892   corrMB->SetLineColor(2);
1893   corrMB->SetMarkerStyle(25);
1894   corrMB->SetMarkerColor(2);
1895   corrMB->Draw("SAME PE");
1896   
1897   corrNSD->SetLineColor(4);
1898   corrNSD->SetMarkerStyle(24);
1899   corrNSD->SetMarkerColor(4);
1900   corrNSD->Draw("SAME PE");
1901   
1902   Printf("       MB  INEL  NSD");
1903   Printf("bin 0: %f %f %f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1));
1904   Printf("bin 1: %f %f %f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2));
1905
1906   TLegend* legend = new TLegend(0.3, 0.6, 0.85, 0.85);
1907   legend->SetFillColor(0);
1908   legend->AddEntry(corrINEL, "Correction to inelastic sample");
1909   legend->AddEntry(corrNSD, "Correction to NSD sample");
1910   legend->AddEntry(corrMB, "Correction to triggered sample");
1911   legend->SetTextSize(0.04);
1912
1913   legend->Draw();
1914
1915   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1916 }
1917
1918 void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1919 {
1920   loadlibs();
1921
1922   TFile::Open(correctionFile);
1923   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1924   mult->LoadHistograms("Multiplicity");
1925
1926   TFile::Open(measuredFile);
1927   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1928   mult2->LoadHistograms("Multiplicity");
1929
1930   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1931
1932   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1933
1934   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 1e5);
1935   
1936   TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1937   
1938   return;
1939   
1940   TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1941
1942   TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1943
1944   if (!mc)
1945   {
1946     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1947     DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1948   }
1949
1950   TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
1951   errorResponse->Write();
1952   errorMeasured->Write();
1953   errorBoth->Write();
1954   file->Close();
1955 }
1956
1957 void DrawStatisticalUncertainty()
1958 {
1959   TFile::Open("StatisticalUncertainty.root");
1960   
1961   errorResponse = (TH1*) gFile->Get("errorResponse");
1962   errorMeasured = (TH1*) gFile->Get("errorMeasured");
1963   errorBoth = (TH1*) gFile->Get("errorBoth");
1964   
1965   TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1966   canvas->SetGridx();
1967   canvas->SetGridy();
1968   canvas->SetRightMargin(0.05);
1969   canvas->SetTopMargin(0.05);
1970
1971   errorResponse->SetLineColor(1);
1972   errorResponse->GetXaxis()->SetRangeUser(0, longDisplayRange);
1973   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1974   errorResponse->SetStats(kFALSE);
1975   errorResponse->SetTitle(";true multiplicity;Uncertainty");
1976
1977   errorResponse->Draw();
1978
1979   errorMeasured->SetLineColor(2);
1980   errorMeasured->Draw("SAME");
1981
1982   errorBoth->SetLineColor(4);
1983   errorBoth->Draw("SAME");
1984
1985   Printf("Average errorResponse: %f", errorResponse->Integral(2, displayRange) / (displayRange - 1));
1986   Printf("Average errorMeasured: %f", errorMeasured->Integral(2, displayRange) /  (displayRange - 1));
1987   Printf("Average errorBoth: %f", errorBoth->Integral(2, displayRange) /  (displayRange - 1));
1988
1989   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1990 }
1991
1992 void StatisticalUncertaintyCompare(const char* det = "SPD")
1993 {
1994   TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1995   TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1996   TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1997   TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1998
1999   TString str;
2000   str.Form("StatisticalUncertaintyCompare%s", det);
2001
2002   TCanvas* canvas = new TCanvas(str, str, 800, 500);
2003   canvas->SetGridx();
2004   canvas->SetGridy();
2005   canvas->SetRightMargin(0.05);
2006   canvas->SetTopMargin(0.05);
2007   
2008   errorResponse->Scale(1.0 / sqrt(2));
2009   errorMeasured->Scale(1.0 / sqrt(2));
2010   errorBoth->Scale(1.0 / sqrt(2));
2011
2012   errorResponse->SetLineColor(1);
2013   errorResponse->GetXaxis()->SetRangeUser(0, displayRange);
2014   errorResponse->GetYaxis()->SetRangeUser(0, 0.18);
2015   errorResponse->SetStats(kFALSE);
2016   errorResponse->GetYaxis()->SetTitleOffset(1.2);
2017   errorResponse->SetTitle(Form(";%s;#sqrt{2}^{-1} #sigma(unfolded - unfolded_{0}) / unfolded_{0}", GetMultLabel()));
2018
2019   errorResponse->Draw();
2020
2021   errorMeasured->SetLineColor(2);
2022   errorMeasured->Draw("SAME");
2023
2024   errorBoth->SetLineColor(4);
2025   errorBoth->Draw("SAME");
2026
2027   TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
2028   TH1* errorResponse2 = (TH1*) file2->Get("errorResponse");
2029   TH1* errorMeasured2 = (TH1*) file2->Get("errorMeasured");
2030   TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
2031
2032   errorResponse2->Scale(1.0 / sqrt(2));
2033   errorMeasured2->Scale(1.0 / sqrt(2));
2034   errorBoth2->Scale(1.0 / sqrt(2));
2035   
2036   errorResponse2->SetLineStyle(2);
2037   errorResponse2->Draw("SAME");
2038   
2039   errorMeasured2->SetLineColor(2);
2040   errorMeasured2->SetLineStyle(2);
2041   errorMeasured2->Draw("SAME");
2042   
2043   errorBoth2->SetLineColor(4);
2044   errorBoth2->SetLineStyle(2);
2045   errorBoth2->Draw("SAME");
2046
2047   TLegend* legend = new TLegend(0.2, 0.5, 0.8, 0.9);
2048   legend->SetFillColor(0);
2049   legend->SetTextSize(0.04);
2050   legend->AddEntry(errorBoth, "Both (Bayesian unfolding)");
2051   legend->AddEntry(errorMeasured, "Measured (Bayesian unfolding)");
2052   legend->AddEntry(errorResponse, "Response matrix (Bayesian unfolding)");
2053   legend->AddEntry(errorBoth2, "Both (#chi^{2}-minimization)");
2054   legend->AddEntry(errorMeasured2, "Measured (#chi^{2}-minimization)");
2055   legend->AddEntry(errorResponse2, "Response matrix (#chi^{2}-minimization)");
2056   legend->Draw();
2057
2058   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2059 }
2060
2061 void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
2062 {
2063  const char* files[] = { "multiplicityMC_nd.root", "multiplicityMC_sd.root", "multiplicityMC_dd.root", "multiplicityMC_xsection.root" };
2064
2065   loadlibs();
2066
2067   TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 600);
2068   canvas->SetGridx();
2069   canvas->SetGridy();
2070   canvas->SetRightMargin(0.05);
2071   canvas->SetTopMargin(0.05);
2072
2073   AliMultiplicityCorrection* data[4];
2074   TH1* effArray[4];
2075   TH1* effErrorArray[2];
2076
2077   Int_t markers[] = { 24, 25, 26, 5 };
2078   //Int_t markers[] = { 2, 25, 24, 5 };
2079   Int_t colors[] = { 1, 2, 4, 6 };
2080   //Int_t colors[] = { 1, 1, 1, 1 };
2081
2082   //TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
2083   TLegend* legend = new TLegend(0.3, 0.3, 0.9, 0.6);
2084   legend->SetTextSize(0.04);
2085   legend->SetFillColor(0);
2086  
2087   for (Int_t i=0; i<4; ++i)
2088   {
2089     TString name;
2090     name.Form("Multiplicity_%d", i);
2091
2092     TFile::Open(files[i]);
2093     data[i] = new AliMultiplicityCorrection(name, name);
2094
2095     if (i < 3)
2096     {
2097       data[i]->LoadHistograms("Multiplicity");
2098     }
2099     else
2100       data[i]->LoadHistograms("Multiplicity_0");
2101
2102     TH1* eff = 0;
2103     if (eventType == -1)
2104     {
2105       eff = (TH1*) data[i]->GetTriggerEfficiency(etaRange)->Clone(Form("eff_%d", i));
2106     }
2107     else
2108       eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
2109     effArray[i] = eff;
2110
2111     eff->GetXaxis()->SetRangeUser(0, 15);
2112     eff->GetYaxis()->SetRangeUser(0, 1.19);
2113     eff->SetStats(kFALSE);
2114     eff->GetXaxis()->SetTitle(GetMultLabel());
2115     eff->GetYaxis()->SetTitle("Efficiency");
2116     eff->SetTitle("");
2117     eff->SetLineColor(colors[i]);
2118     eff->SetMarkerColor(colors[i]);
2119     eff->SetMarkerStyle(markers[i]);
2120
2121     if (i == 3)
2122     {
2123       // once for INEL, once for NSD
2124       for (AliMultiplicityCorrection::EventType eventType2 = AliMultiplicityCorrection::kINEL; eventType2 <= AliMultiplicityCorrection::kNSD; eventType2++)
2125       {
2126         effDiff = (TH1*) data[i]->GetEfficiency(etaRange, eventType2)->Clone(Form("effDiff_%d", i));
2127         
2128         for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2129           effDiff->SetBinError(bin, 0);
2130   
2131         // loop over cross section combinations
2132         for (Int_t j=1; j<7; ++j)
2133         {
2134           AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
2135           mult->LoadHistograms(Form("Multiplicity_%d", j));
2136   
2137           TH1* eff2 = mult->GetEfficiency(etaRange, eventType2);
2138   
2139           for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2140           {
2141             // TODO we could also do asymmetric errors here
2142             Float_t deviation = TMath::Abs(effDiff->GetBinContent(bin) - eff2->GetBinContent(bin));
2143   
2144             effDiff->SetBinError(bin, TMath::Max(effDiff->GetBinError(bin), (Double_t) deviation));
2145           }
2146         }
2147   
2148         for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2149         {
2150           //if (eventType2 == AliMultiplicityCorrection::kINEL)
2151             //eff->SetBinError(bin, 0);
2152             //eff->SetBinError(bin, effDiff->GetBinError(bin));
2153           if (bin < 20 && effDiff->GetBinContent(bin) > 0)
2154             Printf("Bin %d: Error: %.2f", bin, 100.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2155         }
2156         
2157         if (uncertainty) {
2158                 TH1* effError = (TH1*) effDiff->Clone(Form("effError_%s", (eventType2 == AliMultiplicityCorrection::kINEL) ? "INEL" : "NSD"));
2159           effErrorArray[eventType2 - AliMultiplicityCorrection::kINEL] = effError;
2160                 effError->Reset();
2161         
2162                 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2163                   if (effDiff->GetBinContent(bin) > 0)
2164                     effError->SetBinContent(bin, 1.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2165         
2166                 effError->SetLineColor(1);
2167                 effError->SetMarkerStyle(1);
2168                 
2169           if (eventType2 == AliMultiplicityCorrection::kNSD)
2170             effError->SetLineStyle(2);
2171           
2172           effError->DrawCopy("SAME HIST");
2173         }
2174       }
2175     }
2176
2177     canvas->cd();
2178     if (i == 0)
2179     {
2180       eff->DrawCopy("P");
2181     }
2182     else
2183       eff->DrawCopy("SAME P");
2184
2185     legend->AddEntry(eff, (((i == 0) ? "Non-diffractive" : ((i == 1) ? "Single-diffractive" : ((i == 2) ? "Double-diffractive" : "Pythia combined")))));
2186   }
2187
2188   if (uncertainty)
2189   {
2190     legend->AddEntry(effErrorArray[0], "Relative syst. uncertainty: inelastic");
2191     legend->AddEntry(effErrorArray[1], "Relative syst. uncertainty: NSD");
2192   
2193     file = TFile::Open("uncertainty_xsection.root", "RECREATE");
2194     effErrorArray[0]->Write();
2195     effErrorArray[1]->Write();
2196     file->Close();
2197   }
2198
2199   legend->Draw();
2200
2201   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2202 }
2203
2204 void ModelDependencyPlot()
2205 {
2206   loadlibs();
2207
2208   TFile::Open("multiplicityMC.root");
2209   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2210   mult->LoadHistograms("Multiplicity");
2211   
2212   hist = mult->GetCorrelation(etaRange);
2213   
2214   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
2215   {
2216     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
2217     {
2218       hist->SetBinContent(0, y, z, 0);
2219       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2220     }
2221   }
2222
2223   TH2* proj = (TH2*) hist->Project3D("zy");
2224
2225   TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 1200, 600);
2226
2227   canvas->Divide(2, 1);
2228
2229   canvas->cd(2);
2230   gPad->SetLogy();
2231   gPad->SetGridx();
2232   gPad->SetGridy();
2233   gPad->SetTopMargin(0.05);
2234   gPad->SetRightMargin(0.05);
2235  
2236   Int_t selectedMult = 30;
2237   Int_t yMax = 9e4;
2238
2239   TH1* full = proj->ProjectionX("full");
2240   TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult)); 
2241
2242   full->SetStats(kFALSE);
2243   full->GetXaxis()->SetRangeUser(0, displayRange);
2244   full->GetYaxis()->SetRangeUser(5, yMax);
2245   full->SetTitle(";Multiplicity;Entries");
2246
2247   selected->SetLineColor(0);
2248   selected->SetMarkerColor(2);
2249   selected->SetMarkerStyle(5);
2250
2251   full->Draw();
2252   selected->Draw("SAME P");
2253
2254   TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
2255   legend->SetFillColor(0);
2256   legend->AddEntry(full, "True");
2257   legend->AddEntry(selected, "Measured");
2258   legend->Draw();
2259  
2260   TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
2261   line->SetLineWidth(2);
2262   line->Draw();
2263
2264   canvas->cd(1);
2265   gPad->SetLogy();
2266   gPad->SetGridx();
2267   gPad->SetGridy();
2268   gPad->SetTopMargin(0.05);
2269   gPad->SetRightMargin(0.05);
2270
2271   full = proj->ProjectionY("full2");
2272   selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
2273
2274   full->SetStats(kFALSE);
2275   full->GetXaxis()->SetRangeUser(0, displayRange);
2276   full->GetYaxis()->SetRangeUser(5, yMax);
2277   full->SetTitle(";Multiplicity;Entries");
2278
2279   full->SetLineColor(0);
2280   full->SetMarkerColor(2);
2281   full->SetMarkerStyle(5);
2282
2283   full->Draw("P");
2284   selected->Draw("SAME");
2285
2286   legend->Draw();
2287
2288   line = new TLine(selectedMult, 5, selectedMult, yMax);
2289   line->SetLineWidth(2);
2290   line->Draw();
2291
2292   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2293 }
2294
2295 void SystematicpTSpectrum()
2296 {
2297   gSystem->Load("libPWG0base");
2298
2299   TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
2300   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2301   mult->LoadHistograms("Multiplicity");
2302
2303   TFile::Open("multiplicityMC_100k_syst.root");
2304   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2305   mult2->LoadHistograms("Multiplicity");
2306
2307   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2308   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2309   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2310
2311   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2312   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2313
2314   DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
2315 }
2316
2317 void FitPt(const char* fileName)
2318 {
2319   // needs a MC file from the dNdEta analysis
2320
2321   TFile::Open(fileName);
2322
2323   TH1* genePt = (TH1*) gFile->Get("fdNdpT");
2324   
2325   genePt->SetTitle(";p_{T} (GeV/c);1/p_{T} dN_{ch}/dp_{T} (GeV/c)^{-2}");
2326   // number of events
2327   genePt->Scale(1.0 / 287800);
2328   // bin width
2329   genePt->Scale(1.0 / genePt->GetXaxis()->GetBinWidth(1));
2330   
2331   genePt->GetXaxis()->SetRangeUser(0, 0.4);
2332
2333   TF1* func = new TF1("func", "[1]*x*exp(x*[0])");
2334   func->SetParameters(-1, 1);
2335
2336   genePt->SetMarkerStyle(25);
2337   genePt->SetTitle("");
2338   genePt->SetStats(kFALSE);
2339
2340   new TCanvas;
2341   genePt->DrawCopy("P");
2342   func->DrawCopy("SAME");
2343   gPad->SetLogy();
2344
2345   genePt->Fit(func, "0", "", 0, 0.25);
2346   genePt->Fit(func, "0", "", 0, 0.25);
2347
2348   TCanvas* canvas = new TCanvas("FitPt", "FitPt", 600, 600);
2349
2350   gPad->SetGridx();
2351   gPad->SetGridy();
2352   gPad->SetLeftMargin(0.13);
2353   gPad->SetRightMargin(0.05);
2354   gPad->SetTopMargin(0.05);
2355
2356   //genePt->GetXaxis()->SetRangeUser(0, 1);
2357   genePt->GetYaxis()->SetRangeUser(2, 200);
2358   genePt->GetYaxis()->SetTitleOffset(1.4);
2359   genePt->GetXaxis()->SetTitleOffset(1.1);
2360   genePt->DrawCopy("P");
2361   //func->SetRange(0, 0.3);
2362   func->DrawCopy("SAME");
2363   gPad->SetLogy();
2364
2365   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2366   
2367   TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2368
2369   TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2370
2371   TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2372
2373   for (Int_t param=0; param<2; param++)
2374   {
2375     for (Int_t sign=0; sign<2; sign++)
2376     {
2377       TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign));  // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2378       func2->SetParameters(func->GetParameters());
2379       //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2380
2381       Float_t factor = ((sign == 0) ? 0.75 : 1.25);
2382       func2->SetParameter(param, func2->GetParameter(param) * factor);
2383       //func2->Print();
2384
2385       canvas->cd();
2386       func2->SetLineWidth(2);
2387       func2->SetLineColor(2);
2388       func2->SetLineStyle(2);
2389       func2->DrawCopy("SAME");
2390
2391       canvas2->cd();
2392       TH1* second = func2->GetHistogram();
2393       second->Divide(first);
2394       second->SetLineColor(param + 1);
2395       // set to 1 above 0.2 GeV
2396       for (Int_t bin=second->GetXaxis()->FindBin(0.20001); bin<=second->GetNbinsX(); bin++)
2397         second->SetBinContent(bin, 1);
2398       second->GetYaxis()->SetRangeUser(0, 2);
2399       second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2400       second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2401     }
2402   }
2403
2404   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2405   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2406
2407   file->Close();
2408 }
2409
2410 void DrawSystematicpT()
2411 {
2412   TFile* file = TFile::Open("SystematicpT.root");
2413
2414   TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2415   TH1* result2 = (TH1*) file->Get("result_unity");
2416
2417   TH1* mcHist[12];
2418   TH1* result[12];
2419
2420   Int_t nParams = 5;
2421
2422   for (Int_t id=0; id<nParams*2; ++id)
2423   {
2424     mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2425     result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2426   }
2427
2428   DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2429
2430   //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2431
2432   DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2433
2434   //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2435
2436   // does not make sense: mc is different
2437   //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2438 }
2439
2440 void SystematicpT(Bool_t chi2 = 1)
2441 {
2442   gSystem->Load("libPWG0base");
2443
2444   TFile::Open("ptspectrum900.root");
2445   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2446   mult->LoadHistograms("Multiplicity");
2447
2448   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2449
2450   TH1* mcHist[12];
2451   TH1* result[12];
2452
2453   Int_t nParams = 5;
2454
2455   for (Int_t param=0; param<nParams; param++)
2456   {
2457     for (Int_t sign=0; sign<2; sign++)
2458     {
2459       // calculate result with systematic effect
2460       TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2461       mult2->LoadHistograms("Multiplicity");
2462
2463       mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2464
2465       if (chi2)
2466       {
2467         mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2468         mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2469       }
2470       else
2471         mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2472
2473       Int_t id = param * 2 + sign;
2474
2475       mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2476       result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2477
2478       TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2479       DrawResultRatio(mcHist[id], result[id], tmp);
2480     }
2481   }
2482
2483   // calculate normal result
2484   TFile::Open("ptspectrum100_1.root");
2485   mult2->LoadHistograms("Multiplicity");
2486   TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2487
2488   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2489
2490   if (chi2)
2491   {
2492     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2493   }
2494   else
2495     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2496
2497   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2498
2499   TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2500   mcHist2->Write();
2501   result2->Write();
2502   for (Int_t id=0; id<nParams*2; ++id)
2503   {
2504     mcHist[id]->Write();
2505     result[id]->Write();
2506   }
2507   file->Close();
2508
2509   DrawSystematicpT();
2510 }
2511
2512 void DrawSystematicpT2()
2513 {
2514   //displayRange = 200;
2515
2516   // read from file
2517   TFile* file = TFile::Open("SystematicpT2.root");
2518   TH1* mcHist = (TH1*) file->Get("mymc");
2519   TH1* result[12];
2520   result[0] = (TH1*) file->Get("result_unity");
2521   Int_t nParams = 5;
2522   for (Int_t id=0; id<nParams*2; ++id)
2523     result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2524
2525   DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2526   DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2527   DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2528 }
2529
2530 void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2531 {
2532   gSystem->Load("libPWG0base");
2533
2534   if (tpc)
2535   {
2536     SetTPC();
2537     TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2538   }
2539   else
2540     TFile::Open("ptspectrum100_1.root");
2541
2542   AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2543   measured->LoadHistograms("Multiplicity");
2544   TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2545
2546   TH1* result[12];
2547
2548   Int_t nParams = 5;
2549
2550   // -1 = unity change, 0...4 parameters
2551   for (Int_t id=-1; id<nParams*2; id++)
2552   {
2553     Int_t param = id / 2;
2554     Int_t sign = id % 2;
2555
2556     TString idStr;
2557     if (id == -1)
2558     {
2559       idStr = "unity";
2560     }
2561     else
2562       idStr.Form("%d_%d", param, sign);
2563
2564     // calculate result with systematic effect
2565     if (tpc)
2566     {
2567       TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2568     }
2569     else
2570       TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2571
2572     AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2573     response->LoadHistograms("Multiplicity");
2574
2575     response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2576
2577     if (chi2)
2578     {
2579       response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2580       response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2581     }
2582     else
2583       response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2584
2585     result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2586
2587     TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2588     DrawResultRatio(mcHist, result[id+1], tmp);
2589   }
2590
2591   TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2592   mcHist->Write();
2593   for (Int_t id=0; id<nParams*2+1; ++id)
2594     result[id]->Write();
2595   file->Close();
2596
2597   DrawSystematicpT2();
2598 }
2599
2600 void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2601 {
2602   // only needed for TPC
2603   SetTPC();
2604
2605   gSystem->Load("libPWG0base");
2606
2607   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2608   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2609   mult->LoadHistograms("Multiplicity");
2610
2611   TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2612   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2613   mult2->LoadHistograms("Multiplicity");
2614
2615   // "normal" result
2616   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2617
2618   if (chi2)
2619   {
2620     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2621     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2622   }
2623   else
2624     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2625
2626   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2627   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2628
2629   TH1* syst[2];
2630
2631   // change of pt spectrum (down)
2632   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2633   AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2634   mult3->LoadHistograms("Multiplicity");
2635   mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2636   if (chi2)
2637   {
2638     mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2639     mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2640   }
2641   else
2642     mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2643   syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2644
2645   // change of pt spectrum (up)
2646   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2647   AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2648   mult4->LoadHistograms("Multiplicity");
2649   mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2650   if (chi2)
2651   {
2652     mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2653     mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2654   }
2655   else
2656     mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2657   syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
2658
2659   DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
2660
2661   Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2662   Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2663 }
2664
2665 TH1* SystematicsSummary(Bool_t tpc = 0, Bool_t nsd = kTRUE)
2666 {
2667   Int_t nEffects = 7;
2668
2669   TH1* effects[10];
2670   const char** names = 0;
2671   Int_t colors[] = { 1, 2, 4, 1, 2, 4 };
2672   Int_t styles[] = { 1, 2, 3, 1, 2, 3 };
2673   Int_t widths[] = { 1, 1, 1, 2, 2, 2 };
2674   Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
2675   
2676   TH1* dummy = new TH2F("dummy", Form(";%s;Uncertainty", GetMultLabel()), 202, -1.5, 200.5, 100, 0, 0.4);
2677   dummy->SetStats(0);
2678
2679   for (Int_t i=0; i<nEffects; ++i)
2680     effects[i] = new TH1F("SystematicsSummary", Form(";%s;Uncertainty", GetMultLabel()), 201, -0.5, 200.5);
2681
2682   if (tpc)
2683   {
2684     SetTPC();
2685
2686     const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2687     names = namesTPC;
2688
2689     // method
2690     TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2691     TH1* hist = (TH1*) file->Get("errorBoth");
2692
2693     // smooth a bit, but skip 0 bin
2694     effects[0]->SetBinContent(2, hist->GetBinContent(2));
2695     for (Int_t i=3; i<=200; ++i)
2696       effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2697
2698     // relative x-section
2699     effects[1]->SetBinContent(2, 0.005);
2700     effects[1]->SetBinContent(3, 0.0025);
2701     effects[1]->SetBinContent(4, 0.0025);
2702
2703     // particle composition
2704     for (Int_t i=2; i<=101; ++i)
2705     {
2706       if (i < 41)
2707       {
2708         effects[2]->SetBinContent(i, 0.01);
2709       }
2710       else if (i < 76)
2711       {
2712         effects[2]->SetBinContent(i, 0.02);
2713       }
2714       else
2715         effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2716     }
2717
2718     // pt cut off (only tpc)
2719     for (Int_t i=2; i<=101; ++i)
2720     {
2721       if (i < 11)
2722       {
2723         effects[3]->SetBinContent(i, 0.05);
2724       }
2725       else if (i < 51)
2726       {
2727         effects[3]->SetBinContent(i, 0.01);
2728       }
2729       else
2730         effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2731     }
2732
2733     // track selection (only tpc)
2734     for (Int_t i=2; i<=101; ++i)
2735       effects[4]->SetBinContent(i, 0.03);
2736
2737     // secondaries
2738     for (Int_t i=2; i<=101; ++i)
2739       effects[5]->SetBinContent(i, 0.01);
2740
2741     // pt spectrum
2742     for (Int_t i=2; i<=101; ++i)
2743     {
2744       if (i < 21)
2745       {
2746         effects[6]->SetBinContent(i, 0.05);
2747       }
2748       else if (i < 51)
2749       {
2750         effects[6]->SetBinContent(i, 0.02);
2751       }
2752       else
2753         effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2754     }
2755
2756   }
2757   else
2758   {
2759     nEffects = 5;
2760
2761     //const char* namesSPD[] = { "Particle composition",  "p_{t} cut-off", "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)"};
2762     const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)", "Particle composition",  "p_{t} cut-off"};
2763     names = namesSPD;
2764
2765     currentEffect = 0;
2766
2767     // method
2768     TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2769     TH1* hist = (TH1*) file->Get("errorBoth");
2770     hist->Scale(1.0 / sqrt(2));
2771
2772     // smooth a bit, but skip 0 bin
2773     /*effects[currentEffect]->SetBinContent(1, hist->GetBinContent(1));
2774     for (Int_t i=2; i<=201; ++i)
2775       effects[currentEffect]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);*/
2776     effects[currentEffect] = hist;
2777
2778     currentEffect++;
2779
2780     // relative x-section
2781     file = TFile::Open("uncertainty_xsection.root");
2782     effects[currentEffect++] = (TH1*) file->Get("effError_INEL");
2783     effects[currentEffect] = (TH1*) file->Get("effError_NSD");
2784     effects[currentEffect]->SetLineStyle(1);
2785     //effects[2]->SetBinContent(1, 0.20);
2786     //effects[2]->SetBinContent(2, 0.01);
2787     //effects[2]->SetBinContent(3, 0.002);
2788     
2789     currentEffect++;
2790     
2791     // particle composition
2792     effects[currentEffect]->SetBinContent(1, 0.16);
2793     for (Int_t i=2; i<=81; ++i)
2794     {
2795       effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * i / 81);
2796     }
2797     
2798     currentEffect++;
2799
2800     // pt spectrum
2801     effects[currentEffect]->SetBinContent(1, 0.06);
2802     effects[currentEffect]->SetBinContent(2, 0.03);
2803     for (Int_t i=3; i<=81; ++i)
2804     {
2805       if (i <= 61)
2806       {
2807         effects[currentEffect]->SetBinContent(i, 0.01);
2808       }
2809       else if (i <= 81)
2810       {
2811         effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * (i - 61) / 20);
2812       }
2813     }
2814     
2815 //     currentEffect++;
2816 //         
2817 //     // material budget
2818 //     for (Int_t i=1; i<=81; ++i)
2819 //     {
2820 //       if (i < 5)
2821 //         effects[currentEffect]->SetBinContent(i, 0.05 - 0.01 * i);
2822 //       if (i > 51)
2823 //         effects[currentEffect]->SetBinContent(i, 0.05 * (i - 50) / 30);
2824 //     }
2825 //         
2826     currentEffect++;
2827     
2828   }
2829
2830   TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 500);
2831   canvas->SetRightMargin(0.05);
2832   canvas->SetTopMargin(0.05);
2833   //canvas->SetGridx();
2834   canvas->SetGridy();
2835   TLegend* legend = new TLegend(0.2, 0.4, 0.7, 0.4 + 0.5 * nEffects / 7);
2836   legend->SetFillColor(0);
2837   legend->SetTextSize(0.04);
2838   dummy->Draw();
2839   dummy->GetXaxis()->SetRangeUser(0, displayRange);
2840
2841   for (Int_t i=0; i<nEffects; ++i)
2842   {
2843     TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2844     /*current->Reset();
2845     for (Int_t j=0; j<nEffects-i; ++j)
2846       current->Add(effects[j]);*/
2847
2848     current->SetLineColor(colors[i]);
2849     current->SetLineStyle(styles[i]);
2850     current->SetLineWidth(widths[i]);
2851     //current->SetFillColor(colors[i]);
2852     current->SetMarkerColor(colors[i]);
2853     //current->SetMarkerStyle(markers[i]);
2854
2855     current->SetStats(kFALSE);
2856     current->GetYaxis()->SetRangeUser(0, 0.4);
2857     current->DrawCopy("SAME");
2858     legend->AddEntry(current, names[i]);
2859
2860     //TLatex* text = new TLatex(displayRange+2, current->GetBinContent(displayRange+1), names[i]);
2861     TLatex* text = new TLatex(displayRange+2, 0.1 - i * 0.02, names[i]);
2862     text->SetTextSize(0.04);
2863     text->SetTextColor(colors[i]);
2864     //text->Draw();
2865   }
2866
2867   // add total in square
2868   TH1* totalINEL = (TH1*) effects[0]->Clone("totalINEL");
2869   totalINEL->Reset();
2870   TH1* totalNSD = (TH1*) totalINEL->Clone("totalNSD");
2871
2872   for (Int_t i=0; i<nEffects; ++i)
2873   {
2874     //Printf("%d %f", i, effects[i]->GetBinContent(20));
2875     effects[i]->Multiply(effects[i]);
2876     
2877     if (i != 2)
2878       totalINEL->Add(effects[i]);
2879     if (i != 1)
2880       totalNSD->Add(effects[i]);
2881   }
2882   
2883   for (Int_t i=1; i<=totalINEL->GetNbinsX(); ++i)
2884   {
2885     totalINEL->SetBinError(i, 0);
2886     if (totalINEL->GetBinContent(i) > 0)
2887       totalINEL->SetBinContent(i, TMath::Min(sqrt(totalINEL->GetBinContent(i)), 1.0));
2888     totalNSD->SetBinError(i, 0);
2889     if (totalNSD->GetBinContent(i) > 0)
2890       totalNSD->SetBinContent(i, TMath::Min(sqrt(totalNSD->GetBinContent(i)), 1.0));
2891   }
2892
2893   //Printf("%f", total->GetBinContent(20));
2894
2895   totalINEL->SetMarkerStyle(5);
2896   totalINEL->SetMarkerColor(1);
2897   legend->AddEntry(totalINEL, "Total (INEL)", "P");
2898   
2899   totalNSD->SetMarkerStyle(24);
2900   totalNSD->SetMarkerColor(2);
2901   legend->AddEntry(totalNSD, "Total (NSD)", "P");
2902   
2903   Printf("total in bin 0 is INEL: %f NSD: %f", totalINEL->GetBinContent(1), totalNSD->GetBinContent(1));
2904   totalINEL->DrawCopy("SAME P"); //->SetBinContent(1, 0);
2905   totalNSD->DrawCopy("SAME P"); //->SetBinContent(1, 0);
2906
2907   legend->Draw();
2908
2909   canvas->SaveAs(canvas->GetName());
2910   
2911   file = TFile::Open(Form("%s_syst_error.root", (tpc) ? "tpc" : "spd"), "RECREATE");
2912   totalINEL->Write("inel_1");
2913   totalNSD->Write("nsd_1");
2914   file->Close();
2915
2916   return (nsd) ? totalNSD : totalINEL;
2917 }
2918
2919 void finalPlot(Bool_t tpc = 0, Bool_t small = kFALSE)
2920 {
2921   loadlibs();
2922
2923   if (tpc)
2924     SetTPC();
2925
2926   //TH1* errorNSD = SystematicsSummary(tpc, 1);
2927
2928   TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", 800, 500); 
2929   canvas->SetRightMargin(0.05);
2930   canvas->SetTopMargin(0.05);
2931   canvas->SetGridx();
2932   canvas->SetGridy();
2933   
2934   legend = new TLegend(0.5, 0.6, 0.9, 0.8);
2935   legend->SetFillColor(0);
2936   legend->SetTextSize(0.04);
2937   
2938   for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
2939   {
2940     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
2941     TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
2942     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2943   
2944     DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
2945
2946     // normalize result
2947     //result->Scale(1.0 / result->Integral(2, displayRange));
2948   
2949     result->GetXaxis()->SetRangeUser(0, displayRange);
2950     //result->SetBinContent(1, 0); result->SetBinError(1, 0);
2951     result->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (etaRange+1) * 0.5));
2952     result->SetMarkerStyle(0);
2953     result->SetLineColor(1);
2954     result->SetStats(kFALSE);
2955   
2956     // systematic error
2957     TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
2958     
2959     TH1* systError = (TH1*) result->Clone("systError");
2960     for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
2961       systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2962   
2963     // change error drawing style
2964     systError->SetFillColor(15);
2965     
2966     if (eventType == AliMultiplicityCorrection::kNSD)
2967     {
2968       result->SetLineColor(2);
2969       result->SetMarkerColor(2);
2970       result->SetMarkerStyle(5);
2971     }
2972     
2973     canvas->cd();
2974     systError->DrawCopy(Form("E2 ][ %s", (eventType == AliMultiplicityCorrection::kINEL) ? "" : "SAME"));
2975     result->DrawCopy("SAME E ][");
2976     canvas->SetLogy();
2977     
2978     legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic cross-section" : "NSD cross-section", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
2979   }
2980   
2981   legend->Draw();
2982   /*
2983   //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
2984   TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
2985   text->SetFillColor(0);
2986   text->SetTextAlign(12);
2987   text->AddText("Systematic errors summed quadratically");
2988   text->AddText("0.6 million minimum bias events");
2989   text->AddText("corrected to inelastic events");
2990   text->Draw("B");
2991
2992   TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
2993   text2->SetFillColor(0);
2994   text2->SetTextAlign(12);
2995   text2->AddText("#sqrt{s} = 14 TeV");
2996   if (tpc)
2997   {
2998     text2->AddText("|#eta| < 0.9");
2999   }
3000   else
3001     text2->AddText("|#eta| < 2.0");
3002   text2->AddText("simulated data (PYTHIA)");
3003   text2->Draw("B");
3004   
3005   if (tpc)
3006   {
3007     TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
3008     text3->SetNDC();
3009     text3->Draw();
3010   }
3011   else
3012   {
3013     TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
3014     text3->SetNDC();
3015     text3->Draw();
3016   }
3017
3018   // alice logo
3019   TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
3020   pad->Draw();
3021   pad->cd();
3022   TImage* img = TImage::Open("$HOME/alice.png");
3023   img->SetImageQuality(TAttImage::kImgBest);
3024   img->Draw();
3025
3026   canvas->Modified();
3027   */
3028
3029 /*  TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
3030   text->SetTextSize(0.04);
3031   text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
3032   text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
3033   text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
3034   text->Draw();*/
3035
3036
3037   canvas->SaveAs(canvas->GetName());
3038 }
3039
3040 void finalPlot2(Bool_t tpc = 0)
3041 {
3042   loadlibs();
3043
3044   if (tpc)
3045     SetTPC();
3046
3047   //TH1* errorNSD = SystematicsSummary(tpc, 1);
3048
3049   TCanvas* canvas = new TCanvas("finalPlot2.eps", "finalPlot2.eps", 800, 600); 
3050   canvas->SetRightMargin(0.05);
3051   canvas->SetTopMargin(0.05);
3052   canvas->SetGridx();
3053   canvas->SetGridy();
3054   canvas->SetLogy();
3055   
3056   legend = new TLegend(0.5, 0.7, 0.9, 0.9);
3057   legend->SetFillColor(0);
3058   legend->SetTextSize(0.04);
3059   
3060   Int_t displayRanges[] = { 50, 80, 120 };
3061   
3062   TH2* dummy = new TH2F("dummy", ";True multiplicity N_{ch};P(N_{ch})", 100, -0.5, displayRanges[2]+10, 1000, 5e-5, 5);
3063   dummy->SetStats(0);
3064   dummy->Draw();
3065   
3066   TList list;
3067   
3068   // get syst error
3069   file = TFile::Open(Form("%s_syst_error.root", (tpc) ? "tpc" : "spd"));
3070   TH1* totalINEL = (TH1*) file->Get("inel_1");
3071   TH1* totalNSD = (TH1*) file->Get("nsd_1");
3072   
3073   Int_t count = 0;
3074   for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
3075   {
3076     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
3077     //TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
3078     for (Int_t etaR = 0; etaR <= 2; etaR++)
3079     {
3080       TH1* result = mult->GetMultiplicityESDCorrected(etaR);
3081     
3082       //DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
3083   
3084       // normalize result
3085       result->Scale(TMath::Power(5, etaR) / result->Integral(1, displayRanges[etaR]));
3086     
3087       result->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
3088       //result->SetBinContent(1, 0); result->SetBinError(1, 0);
3089       //result->SetTitle(Form(""));
3090       result->SetMarkerStyle(0);
3091       result->SetLineColor(1);
3092       result->SetLineWidth(2);
3093       //result->SetMarkerStyle(4);
3094       //result->SetStats(kFALSE);
3095     
3096       // systematic error
3097       //TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
3098       TH1* error = (eventType == AliMultiplicityCorrection::kNSD) ? totalNSD : totalINEL;
3099       
3100       TH1* systError = (TH1*) result->Clone("systError");
3101       // small hack until we have syst. errors for all eta ranges
3102       Float_t factor = 80.0 / displayRanges[etaR];
3103       for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
3104       {
3105         systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(factor * i));
3106       }
3107     
3108       // change error drawing style
3109       systError->SetFillColor(15);
3110       
3111       if (eventType == AliMultiplicityCorrection::kNSD)
3112       {
3113         result->SetLineColor(2);
3114         result->SetMarkerColor(2);
3115         result->SetMarkerStyle(5);
3116       }
3117       
3118       canvas->cd();
3119       systError->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
3120       systError->DrawCopy("E2 ][ SAME");
3121       list.Add(result);
3122       
3123       if (eventType == AliMultiplicityCorrection::kINEL)
3124       {
3125         TLatex* Tl = new TLatex;
3126         Tl->SetTextSize(0.04);
3127         //Tl->SetBit(TLatex::kTextNDC);
3128         Tl->SetTextAlign(32);
3129
3130         Float_t etaRangeArr[] = { 0.5, 1.0, 1.4 };
3131         TString tmpStr;
3132         tmpStr.Form("|#eta| < %.1f (x %d)", etaRangeArr[etaR], (Int_t) TMath::Power(5, etaR));
3133         if (etaR == 0)
3134           Tl->DrawLatex(36, result->GetBinContent(41), tmpStr);
3135         if (etaR == 1)
3136         {
3137           Tl->SetTextAlign(12);
3138           Tl->DrawLatex(82, result->GetBinContent(81), tmpStr);
3139         }
3140         if (etaR == 2)
3141         {
3142           Tl->SetTextAlign(12);
3143           Tl->DrawLatex(106, result->GetBinContent(101), tmpStr);
3144         }
3145       }
3146       
3147       if (etaR == 0)
3148         legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic events" : "NSD events", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
3149       
3150       count++;
3151     }
3152   }
3153   
3154   for (Int_t i=0; i<list.GetEntries(); i++)
3155     ((TH1*) list.At(i))->DrawCopy("SAME E ][");
3156   
3157   legend->Draw();
3158
3159   canvas->SaveAs(canvas->GetName());
3160 }
3161
3162 TMatrixD* NonInvertable()
3163 {
3164   const Int_t kSize = 5;
3165
3166   TMatrixD matrix(kSize, kSize);
3167   for (Int_t x=0; x<kSize; x++)
3168   {
3169     for (Int_t y=0; y<kSize; y++)
3170     {
3171       if (x == y)
3172       {
3173         if (x == 0 || x == kSize -1)
3174         {
3175           matrix(x, y) = 0.75;
3176         }
3177         else
3178           matrix(x, y) = 0.5;
3179       }
3180       else if (TMath::Abs(x - y) == 1)
3181       {
3182         matrix(x, y) = 0.25;
3183       }
3184     }
3185   }
3186
3187   matrix.Print();
3188
3189   //TMatrixD inverted(matrix);
3190   //inverted.Invert();
3191   
3192   //inverted.Print();
3193   
3194   return new TMatrixD(matrix);
3195 }
3196
3197 void BlobelUnfoldingExample()
3198 {
3199   const Int_t kSize = 20;
3200
3201   TMatrixD matrix(kSize, kSize);
3202   for (Int_t x=0; x<kSize; x++)
3203   {
3204     for (Int_t y=0; y<kSize; y++)
3205     {
3206       if (x == y)
3207       {
3208         if (x == 0 || x == kSize -1)
3209         {
3210           matrix(x, y) = 0.75;
3211         }
3212         else
3213           matrix(x, y) = 0.5;
3214       }
3215       else if (TMath::Abs(x - y) == 1)
3216       {
3217         matrix(x, y) = 0.25;
3218       }
3219     }
3220   }
3221
3222   //matrix.Print();
3223
3224   TMatrixD inverted(matrix);
3225   inverted.Invert();
3226
3227   //inverted.Print();
3228
3229   TH1F* inputDist = new TH1F("inputDist", ";t;Entries", kSize, -0.5, (Float_t) kSize - 0.5);
3230   TVectorD inputDistVector(kSize);
3231   TH1F* unfolded = inputDist->Clone("unfolded");
3232   TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
3233   measuredIdealDist->SetTitle(";m;Entries");
3234   TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
3235
3236   TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
3237   // norm: 1/(sqrt(2pi)sigma)
3238   gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
3239   //gaus->Print();
3240
3241   for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
3242   {
3243     Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
3244     inputDist->SetBinContent(x, value);
3245     inputDistVector(x-1) = value;
3246   }
3247
3248   TVectorD measuredDistIdealVector = matrix * inputDistVector;
3249   
3250   for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
3251     measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
3252
3253   measuredDist->FillRandom(measuredIdealDist, 10000);
3254
3255   // fill error matrix before scaling
3256   TMatrixD covarianceMatrixMeasured(kSize, kSize);
3257   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3258     covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
3259
3260   TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
3261   //covarianceMatrix.Print();
3262
3263   TVectorD measuredDistVector(kSize);
3264   for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
3265     measuredDistVector(x-1) = measuredDist->GetBinContent(x);
3266
3267   TVectorD unfoldedVector = inverted * measuredDistVector;
3268   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3269     unfolded->SetBinContent(x, unfoldedVector(x-1));
3270
3271   TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1200, 600);
3272   canvas->SetTopMargin(0.05);
3273   canvas->Divide(2, 1);
3274
3275   canvas->cd(1);
3276   canvas->cd(1)->SetLeftMargin(0.15);
3277   canvas->cd(1)->SetRightMargin(0.05);
3278   canvas->cd(1)->SetTopMargin(0.05);
3279   gPad->SetGridx();
3280   gPad->SetGridy();
3281   measuredDist->GetYaxis()->SetRangeUser(-600, 2799);
3282   measuredDist->GetYaxis()->SetTitleOffset(1.9);
3283   measuredDist->SetStats(0);
3284   measuredDist->DrawCopy();
3285   gaus->Draw("SAME");
3286
3287   canvas->cd(2);
3288   canvas->cd(2)->SetLeftMargin(0.15);
3289   canvas->cd(2)->SetRightMargin(0.05);
3290   canvas->cd(2)->SetTopMargin(0.05);
3291   gPad->SetGridx();
3292   gPad->SetGridy();
3293   unfolded->GetYaxis()->SetRangeUser(-600, 2799);
3294   unfolded->GetYaxis()->SetTitleOffset(1.9);
3295   unfolded->SetStats(0);
3296   unfolded->DrawCopy();
3297   gaus->Draw("SAME");
3298
3299   canvas->SaveAs("BlobelUnfoldingExample.eps");
3300   
3301   return;
3302   
3303   // now unfold this with Bayesian
3304   loadlibs();
3305   
3306   // fill a multiplicity object
3307   mult = new AliMultiplicityCorrection("mult", "mult");
3308   for (Int_t x=0; x<kSize; x++)
3309   {
3310     mult->GetMultiplicityVtx(0)->SetBinContent(1, x+1, inputDistVector(x));
3311     mult->GetMultiplicityESD(0)->SetBinContent(1, x+1, measuredDistVector(x)*10000);
3312     for (Int_t y=0; y<kSize; y++)
3313       mult->GetCorrelation(0)->SetBinContent(1, x+1, y+1, matrix(x, y));
3314   }
3315   
3316   //mult->DrawHistograms();
3317   
3318   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 0);
3319   //mult->SetCreateBigBin(kFALSE);
3320   mult->ApplyMinuitFit(0, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); //hist2->ProjectionY("mymchist"));
3321   
3322   //mult->ApplyBayesianMethod(0, kFALSE, AliMultiplicityCorrection::kTrVtx, 0, -1, 0, kFALSE);
3323   
3324   mult->DrawComparison("BlobelExample", 0, kFALSE, kTRUE, mult->GetMultiplicityVtx(0)->ProjectionY("mcmchist", 1, mult->GetMultiplicityVtx(0)->GetNbinsX()));
3325   
3326   
3327 }
3328
3329 void E735Fit()
3330 {
3331   TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
3332   fCurrentESD->Sumw2();
3333
3334   // Open the input stream
3335   ifstream in;
3336   in.open("e735data.txt");
3337
3338   while(in.good())
3339   {
3340     Float_t x, y, ye;
3341     in >> x >> y >> ye;
3342
3343     //Printf("%f %f %f", x, y, ye);
3344     fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
3345     fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
3346   }
3347
3348   in.close();
3349
3350   //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
3351
3352   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
3353
3354   TF1* 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])");
3355   func->SetParNames("scaling", "averagen", "k");
3356   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
3357   func->SetParLimits(1, 0.001, 1000);
3358   func->SetParLimits(2, 0.001, 1000);
3359   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
3360
3361   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
3362   lognormal->SetParNames("scaling", "mean", "sigma");
3363   lognormal->SetParameters(1, 1, 1);
3364   lognormal->SetParLimits(0, 0, 10);
3365   lognormal->SetParLimits(1, 0, 100);
3366   lognormal->SetParLimits(2, 1e-3, 10);
3367
3368   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
3369   fCurrentESD->SetStats(kFALSE);
3370   fCurrentESD->SetMarkerStyle(0);
3371   fCurrentESD->SetLineColor(1);
3372   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
3373   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
3374   fCurrentESD->Draw("");
3375   fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
3376   fCurrentESD->Fit(func, "0", "", 0, 150);
3377   func->SetRange(0, 250);
3378   func->Draw("SAME");
3379   printf("chi2 = %f\n", func->GetChisquare());
3380
3381   fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
3382   lognormal->SetLineColor(2);
3383   lognormal->SetLineStyle(2);
3384   lognormal->SetRange(0, 250);
3385   lognormal->Draw("SAME");
3386
3387   gPad->SetLogy();
3388
3389   canvas->SaveAs("E735Fit.eps");
3390 }
3391
3392 void DifferentSamples()
3393 {
3394   loadlibs();
3395
3396   Int_t n = 2;
3397   const char* filesChi2[] = { "chi2_100k_1.root", "chi2_100k_2.root" };
3398   const char* filesBayesian[] = { "bayesian_100k_1.root", "bayesian_100k_2.root" };
3399
3400   TCanvas* canvas = new TCanvas("DifferentSamples", "DifferentSamples", 1200, 600);
3401   canvas->Divide(2, 1);
3402   
3403   legend = new TLegend(0.15, 0.7, 0.65, 0.9);
3404   legend->SetFillColor(0);
3405   legend->SetTextSize(0.04);
3406
3407   for (Int_t i=0; i<n; i++)
3408   {
3409     AliMultiplicityCorrection* chi2 = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3410     TFile::Open(filesChi2[i]);
3411     chi2->LoadHistograms("Multiplicity");
3412
3413     AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3414     TFile::Open(filesBayesian[i]);
3415     bayesian->LoadHistograms("Multiplicity");
3416     
3417     chi2Hist = chi2->GetMultiplicityESDCorrected(etaRange);
3418     bayesianHist = bayesian->GetMultiplicityESDCorrected(etaRange);
3419     
3420     mc = chi2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
3421     
3422     // normalize and divide
3423     chi2Hist->Scale(1.0 / chi2Hist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3424     bayesianHist->Scale(1.0 / bayesianHist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3425     
3426     chi2Hist->Divide(mc, chi2Hist);
3427     bayesianHist->Divide(mc, bayesianHist);
3428     
3429     canvas->cd(i+1);
3430     gPad->SetTopMargin(0.05);
3431     gPad->SetRightMargin(0.05);
3432     //gPad->SetLeftMargin(0.12);
3433     gPad->SetGridx();
3434     gPad->SetGridy(); 
3435     
3436     chi2Hist->GetXaxis()->SetRangeUser(0, displayRange);
3437     chi2Hist->GetYaxis()->SetTitleOffset(1.3);
3438     chi2Hist->SetStats(0);
3439     chi2Hist->SetTitle(Form(";%s;MC / unfolded", GetMultLabel()));
3440     chi2Hist->GetYaxis()->SetRangeUser(0.2, 1.8);
3441     chi2Hist->Draw("HIST");
3442     
3443     for (Int_t x=1; x<=bayesianHist->GetNbinsX(); x++)
3444       bayesianHist->SetBinError(x, 1e-6);
3445     
3446     bayesianHist->SetLineColor(2);
3447     bayesianHist->SetMarkerColor(2);
3448     bayesianHist->SetMarkerStyle(5);
3449     bayesianHist->Draw("HIST E SAME");
3450     
3451     if (i == 0)
3452     {
3453       legend->AddEntry(chi2Hist, "#chi^{2}-minimization", "L");
3454       legend->AddEntry(bayesianHist, "Bayesian unfolding", "LP");
3455     }
3456     legend->Draw();
3457   }
3458   
3459   canvas->SaveAs("DifferentSamples.eps");
3460 }
3461
3462 void PileUp()
3463 {
3464   loadlibs();
3465   
3466   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
3467   hist2d = mult->GetMultiplicityMC(etaRange, AliMultiplicityCorrection::kINEL);
3468   mult1 = hist2d->ProjectionY("mult1", 1, hist2d->GetNbinsX());
3469   
3470   conv = (TH1*) mult1->Clone("conv");
3471   conv->Reset();
3472   
3473   mult1->Scale(1.0 / mult1->Integral());
3474   
3475   for (Int_t i=1; i<=mult1->GetNbinsX(); i++)
3476     for (Int_t j=1; j<=mult1->GetNbinsX(); j++)
3477       conv->Fill(mult1->GetBinCenter(i)+mult1->GetBinCenter(j), mult1->GetBinContent(i) * mult1->GetBinContent(j));
3478   
3479   conv->Scale(1.0 / conv->Integral());
3480   
3481   c = new TCanvas("c", "c", 800, 500);
3482   gPad->SetLogy();
3483   gPad->SetTopMargin(0.05);
3484   gPad->SetRightMargin(0.05);
3485   mult1->SetTitle(Form(";%s;Probability", GetMultLabel()));
3486   mult1->SetStats(0);
3487   gPad->SetGridx();
3488   gPad->SetGridy();
3489   mult1->Draw();
3490   mult1->GetYaxis()->SetRangeUser(1e-7, 2 * mult1->GetMaximum());
3491   mult1->GetXaxis()->SetRangeUser(0, displayRange);
3492   mult1->GetXaxis()->SetTitleOffset(1.15);
3493   conv->SetLineColor(2);
3494   conv->SetMarkerColor(2);
3495   conv->SetMarkerStyle(5);
3496   conv->DrawCopy("SAME P");
3497   
3498   conv->Scale(0.00058);
3499   conv->DrawCopy("SAME P");
3500   
3501   legend = new TLegend(0.73, 0.73, 0.93, 0.93);
3502   legend->SetFillColor(0);
3503   legend->SetTextSize(0.04);
3504   legend->AddEntry(mult1, "1 collision");
3505   legend->AddEntry(conv, "2 collisions", "P");
3506   legend->Draw();
3507   
3508   c->SaveAs("pileup.eps");
3509
3510   new TCanvas;
3511   conv->Divide(mult1);
3512   conv->Draw();
3513 }
3514
3515 void TestErrorDetermination(Int_t nRandomizations)
3516 {
3517   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, 100);
3518   func->SetParNames("scaling", "averagen", "k");
3519   func->SetParameters(1, 15, 2);
3520   
3521   TF1* func2 = new TF1("nbd2", "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, 100);
3522   func2->SetParNames("scaling", "averagen", "k");
3523   func2->SetParLimits(0, 0.5, 2);
3524   func2->SetParLimits(1, 1, 50);
3525   func2->SetParLimits(2, 1, 10);
3526   func2->SetParameters(1, 15, 2);
3527   //func2->FixParameter(0, 1);
3528   
3529   //new TCanvas; func->Draw("L");
3530   
3531   hist1 = new TH1F("hist1", "", 100, 0.5, 100.5);
3532   hist2 = new TH1F("hist2", "", 100, 0.5, 100.5);
3533   hist1->Sumw2();
3534   
3535   TH1* params[3];
3536   params[0] = new TH1F("param_0", Form("param_%d", 0), 100, 0.95, 1.05);
3537   params[1] = new TH1F("param_1", Form("param_%d", 1), 100, 14, 16);
3538   params[2] = new TH1F("param_2", Form("param_%d", 2), 100, 1.8, 2.2);
3539   
3540   const Int_t nTries = 1000;
3541   for (Int_t i=0; i<nTries; i++)
3542   {
3543     hist1->Reset();
3544     
3545     if (nRandomizations == 1)
3546     {
3547       hist1->FillRandom("nbd", 10000);
3548     }
3549     else if (nRandomizations == 2)
3550     {
3551       hist2->Reset();
3552       hist2->FillRandom("nbd", 10000);
3553       hist1->FillRandom(hist2, 10000);
3554     }
3555     else if (nRandomizations == 3)
3556     {
3557       hist2->Reset();
3558       hist1->FillRandom("nbd", 10000);
3559       hist2->FillRandom(hist1, 10000);
3560       hist1->Reset();
3561       hist1->FillRandom(hist2, 10000);
3562     }
3563     else
3564       return;
3565   
3566     //new TCanvas; hist1->Draw();
3567   
3568     hist1->Scale(1.0 / hist1->Integral());
3569     hist1->Fit(func2, "NQ");
3570     hist1->Fit(func2, "NQ");
3571     for (Int_t j=0; j<3; j++)
3572       params[j]->Fill(func2->GetParameter(j));
3573   }
3574   
3575   for (Int_t j=0; j<3; j++)
3576   {
3577     new TCanvas; params[j]->Draw();
3578     params[j]->Fit("gaus");
3579     Printf("sigma of param %d if %f", j, ((TF1*) params[j]->FindObject("gaus"))->GetParameter(2));
3580   }
3581 }
3582
3583 void DrawRawDistributions(const char* fileName = "multiplicityESD.root")
3584 {
3585   loadlibs();
3586
3587   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileName);
3588   
3589   c = new TCanvas("c", "c", 600, 600);
3590   
3591   dummy = new TH2F("dummy", ";measured multiplicity", 100, -0.5, 149.5, 100, 0.5, 4e4);
3592   dummy->SetStats(0);
3593   dummy->Draw();
3594   gPad->SetGridx();
3595   gPad->SetGridy();
3596   gPad->SetLogy();
3597   
3598   Int_t colors[] = { 1, 2, 4 };
3599   
3600   for (Int_t i=2; i>=0; i--)
3601   {
3602     hist = mult->GetMultiplicityESD(i)->ProjectionY();
3603     
3604     hist->SetLineColor(colors[i]);
3605     hist->DrawCopy("SAME");
3606   }
3607   
3608   
3609 }