439eea72d9981677ee36dd57208c8b8e54faecbd
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / plotsMultiplicity.C
1 //
2 // plots for the note about multplicity measurements
3 //
4
5 #if !defined(__CINT__) || defined(__MAKECINT__)
6
7 #include <TCanvas.h>
8 #include <TPad.h>
9 #include <TH1F.h>
10 #include <TH2F.h>
11 #include <TH3F.h>
12 #include <TLine.h>
13 #include <TF1.h>
14 #include <TSystem.h>
15 #include <TFile.h>
16 #include <TLegend.h>
17 #include <TStopwatch.h>
18 #include <TROOT.h>
19 #include <TGraph.h>
20 #include <TMath.h>
21 #include <TPaveText.h>
22 #include <TImage.h>
23 #include <TLatex.h>
24
25 #include "AliMultiplicityCorrection.h"
26 #include "AliCorrection.h"
27 #include "AliCorrectionMatrix3D.h"
28
29 #endif
30
31 const char* correctionFile = "multiplicityMC_2M.root";
32 const char* measuredFile   = "multiplicityMC_1M_3.root";
33 Int_t etaRange = 3;
34 Int_t displayRange = 200; // axis range
35 Int_t ratioRange = 151;   // range to calculate difference
36 Int_t longDisplayRange = 200;
37
38 const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
39 const char* measuredFileTPC   = "multiplicityMC_TPC_0.6M.root";
40 Int_t etaRangeTPC = 1;
41
42 void SetTPC()
43 {
44   correctionFile = correctionFileTPC;
45   measuredFile = measuredFileTPC;
46   etaRange = etaRangeTPC;
47   displayRange = 100;
48   ratioRange = 76;
49   longDisplayRange = 100;
50 }
51
52 void Smooth(TH1* hist, Int_t windowWidth = 20)
53 {
54   TH1* clone = (TH1*) hist->Clone("clone");
55   for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
56   {
57     Int_t min = TMath::Max(2, bin-windowWidth);
58     Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
59     Float_t average = clone->Integral(min, max) / (max - min + 1);
60
61     hist->SetBinContent(bin, average);
62     hist->SetBinError(bin, 0);
63   }
64
65   delete clone;
66 }
67
68 void responseMatrixPlot()
69 {
70   gSystem->Load("libPWG0base");
71
72   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
73
74   TFile::Open(correctionFile);
75   mult->LoadHistograms("Multiplicity");
76
77   TH1* hist = mult->GetCorrelation(etaRange)->Project3D("zy");
78   hist->SetStats(kFALSE);
79
80   hist->SetTitle(";true multiplicity;measured multiplicity;Entries");
81   hist->GetXaxis()->SetRangeUser(0, longDisplayRange);
82   hist->GetYaxis()->SetRangeUser(0, longDisplayRange);
83
84   TCanvas* canvas = new TCanvas("c1", "c1", 800, 600);
85   canvas->SetRightMargin(0.15);
86   canvas->SetTopMargin(0.05);
87
88   gPad->SetLogz();
89   hist->Draw("COLZ");
90
91   canvas->SaveAs("responsematrix.eps");
92 }
93
94 TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
95 {
96   // normalize unfolded result to mc hist
97   result->Scale(1.0 / result->Integral(2, 200));
98   result->Scale(mcHist->Integral(2, 200));
99
100   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
101   canvas->Range(0, 0, 1, 1);
102
103   TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
104   pad1->Draw();
105
106   TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
107   pad2->Draw();
108
109   pad1->SetRightMargin(0.05);
110   pad2->SetRightMargin(0.05);
111
112   // no border between them
113   pad1->SetBottomMargin(0);
114   pad2->SetTopMargin(0);
115
116   pad1->cd();
117
118   mcHist->GetXaxis()->SetLabelSize(0.06);
119   mcHist->GetYaxis()->SetLabelSize(0.06);
120   mcHist->GetXaxis()->SetTitleSize(0.06);
121   mcHist->GetYaxis()->SetTitleSize(0.06);
122   mcHist->GetYaxis()->SetTitleOffset(0.6);
123
124   mcHist->GetXaxis()->SetRangeUser(0, displayRange);
125
126   mcHist->SetTitle(";true multiplicity;Entries");
127   mcHist->SetStats(kFALSE);
128
129   mcHist->DrawCopy("HIST E");
130   gPad->SetLogy();
131
132   result->SetLineColor(2);
133   result->DrawCopy("SAME HISTE");
134
135   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
136   legend->AddEntry(mcHist, "true distribution");
137   legend->AddEntry(result, "unfolded distribution");
138   legend->SetFillColor(0);
139   legend->Draw();
140
141   pad2->cd();
142   pad2->SetBottomMargin(0.15);
143
144   // calculate ratio
145   mcHist->Sumw2();
146   TH1* ratio = (TH1*) mcHist->Clone("ratio");
147   result->Sumw2();
148   ratio->Divide(ratio, result, 1, 1, "");
149   ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
150   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
151
152   ratio->DrawCopy();
153
154   // get average of ratio
155   Float_t sum = 0;
156   for (Int_t i=2; i<=ratioRange; ++i)
157   {
158     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
159   }
160   sum /= ratioRange-1;
161
162   printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
163
164   TLine* line = new TLine(0, 1, displayRange, 1);
165   line->SetLineWidth(2);
166   line->Draw();
167
168   line = new TLine(0, 1.1, displayRange, 1.1);
169   line->SetLineWidth(2);
170   line->SetLineStyle(2);
171   line->Draw();
172   line = new TLine(0, 0.9, displayRange, 0.9);
173   line->SetLineWidth(2);
174   line->SetLineStyle(2);
175   line->Draw();
176
177   canvas->Modified();
178
179   canvas->SaveAs(epsName);
180
181   return canvas;
182 }
183
184 TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
185 {
186   // draws the 3 plots in the upper plot
187   // draws the ratio between result1 and result2 in the lower plot
188
189   // normalize unfolded result to mc hist
190   result1->Scale(1.0 / result1->Integral(2, 200));
191   result1->Scale(mcHist->Integral(2, 200));
192   result2->Scale(1.0 / result2->Integral(2, 200));
193   result2->Scale(mcHist->Integral(2, 200));
194
195   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
196   canvas->Range(0, 0, 1, 1);
197
198   TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
199   pad1->Draw();
200
201   TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
202   pad2->Draw();
203
204   pad1->SetRightMargin(0.05);
205   pad2->SetRightMargin(0.05);
206
207   // no border between them
208   pad1->SetBottomMargin(0);
209   pad2->SetTopMargin(0);
210
211   pad1->cd();
212
213   mcHist->GetXaxis()->SetLabelSize(0.06);
214   mcHist->GetYaxis()->SetLabelSize(0.06);
215   mcHist->GetXaxis()->SetTitleSize(0.06);
216   mcHist->GetYaxis()->SetTitleSize(0.06);
217   mcHist->GetYaxis()->SetTitleOffset(0.6);
218
219   mcHist->GetXaxis()->SetRangeUser(0, displayRange);
220
221   mcHist->SetTitle(";true multiplicity;Entries");
222   mcHist->SetStats(kFALSE);
223
224   mcHist->DrawCopy("HIST E");
225   gPad->SetLogy();
226
227   result1->SetLineColor(2);
228   result1->DrawCopy("SAME HISTE");
229
230   result2->SetLineColor(4);
231   result2->DrawCopy("SAME HISTE");
232
233   TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
234   legend->AddEntry(mcHist, "true distribution");
235   legend->AddEntry(result1, "unfolded distribution (syst)");
236   legend->AddEntry(result2, "unfolded distribution (normal)");
237   legend->SetFillColor(0);
238   legend->Draw();
239
240   pad2->cd();
241   pad2->SetBottomMargin(0.15);
242
243   result1->GetXaxis()->SetLabelSize(0.06);
244   result1->GetYaxis()->SetLabelSize(0.06);
245   result1->GetXaxis()->SetTitleSize(0.06);
246   result1->GetYaxis()->SetTitleSize(0.06);
247   result1->GetYaxis()->SetTitleOffset(0.6);
248
249   result1->GetXaxis()->SetRangeUser(0, displayRange);
250
251   result1->SetTitle(";true multiplicity;Entries");
252   result1->SetStats(kFALSE);
253
254   // calculate ratio
255   result1->Sumw2();
256   TH1* ratio = (TH1*) result1->Clone("ratio");
257   result2->Sumw2();
258   ratio->Divide(ratio, result2, 1, 1, "");
259   ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
260   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
261
262   ratio->DrawCopy();
263
264   // get average of ratio
265   Float_t sum = 0;
266   for (Int_t i=2; i<=ratioRange; ++i)
267   {
268     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
269   }
270   sum /= ratioRange-1;
271
272   printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
273
274   TLine* line = new TLine(0, 1, displayRange, 1);
275   line->SetLineWidth(2);
276   line->Draw();
277
278   line = new TLine(0, 1.1, displayRange, 1.1);
279   line->SetLineWidth(2);
280   line->SetLineStyle(2);
281   line->Draw();
282   line = new TLine(0, 0.9, displayRange, 0.9);
283   line->SetLineWidth(2);
284   line->SetLineStyle(2);
285   line->Draw();
286
287   canvas->Modified();
288
289   canvas->SaveAs(epsName);
290
291   return canvas;
292 }
293
294 TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE)
295 {
296   // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
297   // a systematic effect
298
299   // normalize results
300   result->Scale(1.0 / result->Integral(2, 200));
301
302   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
303   canvas->SetTopMargin(0.05);
304   canvas->SetRightMargin(0.05);
305
306   result->GetXaxis()->SetRangeUser(0, displayRange);
307   result->GetYaxis()->SetRangeUser(0.55, 1.45);
308   result->SetStats(kFALSE);
309
310   // to get the axis how we want it
311   TH1* dummy = (TH1*) result->Clone("dummy");
312   dummy->Reset();
313   dummy->SetTitle(";true multiplicity;Ratio");
314   dummy->DrawCopy();
315   delete dummy;
316
317   Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
318
319   TLegend* legend = new TLegend(0.2, 0.75, 0.35, 0.95);
320   legend->SetFillColor(0);
321
322   for (Int_t n=0; n<nResultSyst; ++n)
323   {
324     resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
325
326     // calculate ratio
327     TH1* ratio = (TH1*) result->Clone("ratio");
328     ratio->Divide(ratio, resultSyst[n], 1, 1, "");
329     ratio->GetXaxis()->SetRangeUser(1, displayRange);
330
331     if (firstMarker)
332       ratio->SetMarkerStyle(5);
333
334     ratio->SetLineColor(colors[n / 2]);
335     if ((n % 2))
336       ratio->SetLineStyle(2);
337
338     TString drawStr("SAME HIST");
339     if (n == 0 && firstMarker)
340       drawStr = "SAME P";
341     if (errors)
342       drawStr += " E";
343
344     ratio->DrawCopy(drawStr);
345
346     if (legendStrings && legendStrings[n])
347       legend->AddEntry(ratio, legendStrings[n]);
348
349     // get average of ratio
350     Float_t sum = 0;
351     for (Int_t i=2; i<=ratioRange; ++i)
352       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
353     sum /= ratioRange-1;
354
355     printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
356   }
357
358   if (legendStrings)
359     legend->Draw();
360
361   TLine* line = new TLine(0, 1, displayRange, 1);
362   line->SetLineWidth(2);
363   line->Draw();
364
365   line = new TLine(0, 1.1, displayRange, 1.1);
366   line->SetLineWidth(2);
367   line->SetLineStyle(2);
368   line->Draw();
369   line = new TLine(0, 0.9, displayRange, 0.9);
370   line->SetLineWidth(2);
371   line->SetLineStyle(2);
372   line->Draw();
373
374   canvas->SaveAs(epsName);
375   canvas->SaveAs(Form("%s.gif", epsName.Data()));
376
377   return canvas;
378 }
379
380 TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
381 {
382   // draws the ratios of each mc to the corresponding result
383
384   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
385   canvas->SetRightMargin(0.05);
386   canvas->SetTopMargin(0.05);
387
388   for (Int_t n=0; n<nResultSyst; ++n)
389   {
390     // normalize
391     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
392     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
393
394     result[n]->GetXaxis()->SetRangeUser(0, displayRange);
395     result[n]->SetStats(kFALSE);
396
397     // calculate ratio
398     TH1* ratio = (TH1*) result[n]->Clone("ratio");
399     ratio->Divide(mc[n], ratio, 1, 1, "B");
400
401     // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
402     ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
403
404     if (smooth)
405       Smooth(ratio);
406
407     ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
408     ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
409
410     if (dashed)
411     {
412       ratio->SetLineColor((n/2)+1);
413       ratio->SetLineStyle((n%2)+1);
414     }
415     else
416       ratio->SetLineColor(n+1);
417
418     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
419
420     // get average of ratio
421     Float_t sum = 0;
422     for (Int_t i=2; i<=ratioRange; ++i)
423       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
424     sum /= ratioRange-1;
425
426     printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
427   }
428
429   TLine* line = new TLine(0, 1, displayRange, 1);
430   line->SetLineWidth(2);
431   line->Draw();
432
433   line = new TLine(0, 1.1, displayRange, 1.1);
434   line->SetLineWidth(2);
435   line->SetLineStyle(2);
436   line->Draw();
437   line = new TLine(0, 0.9, displayRange, 0.9);
438   line->SetLineWidth(2);
439   line->SetLineStyle(2);
440   line->Draw();
441
442   canvas->Modified();
443
444   canvas->SaveAs(epsName);
445   canvas->SaveAs(Form("%s.gif", epsName.Data()));
446
447   return canvas;
448 }
449
450 TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
451 {
452   // draws the ratios of each mc to the corresponding result
453   // deducts from each ratio the ratio of mcBase / resultBase
454
455   // normalize
456   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
457   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
458
459   // calculate ratio
460   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
461   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
462
463   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
464   canvas->SetRightMargin(0.05);
465   canvas->SetTopMargin(0.05);
466
467   for (Int_t n=0; n<nResultSyst; ++n)
468   {
469     // normalize
470     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
471     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
472
473     result[n]->GetXaxis()->SetRangeUser(0, displayRange);
474     result[n]->SetStats(kFALSE);
475
476     // calculate ratio
477     TH1* ratio = (TH1*) result[n]->Clone("ratio");
478     ratio->Divide(mc[n], ratio, 1, 1, "B");
479     ratio->Add(ratioBase, -1);
480
481     ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
482     ratio->GetYaxis()->SetRangeUser(-1, 1);
483     ratio->SetLineColor(n+1);
484     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
485
486     // get average of ratio
487     Float_t sum = 0;
488     for (Int_t i=2; i<=ratioRange; ++i)
489       sum += TMath::Abs(ratio->GetBinContent(i));
490     sum /= ratioRange-1;
491
492     printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, ratioRange, sum);
493   }
494
495   TLine* line = new TLine(0, 0, displayRange, 0);
496   line->SetLineWidth(2);
497   line->Draw();
498
499   line = new TLine(0, 0.1, displayRange, 0.1);
500   line->SetLineWidth(2);
501   line->SetLineStyle(2);
502   line->Draw();
503   line = new TLine(0, -0.1, displayRange, -0.1);
504   line->SetLineWidth(2);
505   line->SetLineStyle(2);
506   line->Draw();
507
508   canvas->Modified();
509
510   canvas->SaveAs(epsName);
511   canvas->SaveAs(Form("%s.gif", epsName.Data()));
512
513   return canvas;
514 }
515
516 TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
517 {
518   // draws the ratios of each mc to the corresponding result
519   // deducts from each ratio the ratio of mcBase / resultBase
520   // smoothens the ratios by a sliding window
521
522   // normalize
523   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
524   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
525
526   // calculate ratio
527   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
528   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
529
530   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
531   canvas->SetRightMargin(0.05);
532   canvas->SetTopMargin(0.05);
533
534   for (Int_t n=0; n<nResultSyst; ++n)
535   {
536     // normalize
537     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
538     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
539
540     result[n]->GetXaxis()->SetRangeUser(0, displayRange);
541     result[n]->SetStats(kFALSE);
542
543     // calculate ratio
544     TH1* ratio = (TH1*) result[n]->Clone("ratio");
545     ratio->Divide(mc[n], ratio, 1, 1, "B");
546     ratio->Add(ratioBase, -1);
547
548     //new TCanvas; ratio->DrawCopy();
549     // clear 0 bin
550     ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
551
552     Smooth(ratio);
553
554     //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
555
556     canvas->cd();
557     ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
558     ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
559     ratio->SetLineColor((n / 2)+1);
560     ratio->SetLineStyle((n % 2)+1);
561     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
562
563     // get average of ratio
564     Float_t sum = 0;
565     for (Int_t i=2; i<=150; ++i)
566       sum += TMath::Abs(ratio->GetBinContent(i));
567     sum /= 149;
568
569     printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
570   }
571
572   TLine* line = new TLine(0, 0, displayRange, 0);
573   line->SetLineWidth(2);
574   line->Draw();
575
576   line = new TLine(0, 0.1, displayRange, 0.1);
577   line->SetLineWidth(2);
578   line->SetLineStyle(2);
579   line->Draw();
580   line = new TLine(0, -0.1, displayRange, -0.1);
581   line->SetLineWidth(2);
582   line->SetLineStyle(2);
583   line->Draw();
584
585   canvas->Modified();
586
587   canvas->SaveAs(epsName);
588   canvas->SaveAs(Form("%s.gif", epsName.Data()));
589
590   return canvas;
591 }
592
593 void DrawResiduals(TH1* measured, TH1* unfoldedFolded, const char* epsName)
594 {
595   // normalize
596   unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, 200));
597   unfoldedFolded->Scale(measured->Integral(2, 200));
598
599   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
600   canvas->Range(0, 0, 1, 1);
601
602   TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
603   pad1->Draw();
604   pad1->SetGridx();
605   pad1->SetGridy();
606
607   TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
608   pad2->Draw();
609   pad2->SetGridx();
610   pad2->SetGridy();
611
612   TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
613   pad3->SetGridx();
614   pad3->SetGridy();
615   pad3->SetRightMargin(0.05);
616   pad3->SetTopMargin(0.05);
617   pad3->Draw();
618
619   pad1->SetRightMargin(0.05);
620   pad2->SetRightMargin(0.05);
621
622   // no border between them
623   pad1->SetBottomMargin(0);
624   pad2->SetTopMargin(0);
625
626   pad1->cd();
627
628   measured->GetXaxis()->SetLabelSize(0.06);
629   measured->GetYaxis()->SetLabelSize(0.06);
630   measured->GetXaxis()->SetTitleSize(0.06);
631   measured->GetYaxis()->SetTitleSize(0.06);
632   measured->GetYaxis()->SetTitleOffset(0.6);
633
634   measured->GetXaxis()->SetRangeUser(0, 150);
635
636   measured->SetTitle(";measured multiplicity;Entries");
637   measured->SetStats(kFALSE);
638
639   measured->DrawCopy("HIST");
640   gPad->SetLogy();
641
642   unfoldedFolded->SetMarkerStyle(5);
643   unfoldedFolded->SetMarkerColor(2);
644   unfoldedFolded->SetLineColor(0);
645   unfoldedFolded->DrawCopy("SAME P");
646
647   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
648   legend->AddEntry(measured, "measured distribution");
649   legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution");
650   legend->SetFillColor(0);
651   legend->Draw();
652
653   pad2->cd();
654   pad2->SetBottomMargin(0.15);
655
656   // calculate ratio
657   measured->Sumw2();
658   TH1* residual = (TH1*) measured->Clone("residual");
659   unfoldedFolded->Sumw2();
660
661   residual->Add(unfoldedFolded, -1);
662
663   // projection
664   TH1* residualHist = new TH1F("residualHist", ";", 15, -3, 3);
665
666   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
667   {
668     if (measured->GetBinError(i) > 0)
669     {
670       residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
671       residual->SetBinError(i, 1);
672
673       residualHist->Fill(residual->GetBinContent(i));
674     }
675     else
676     {
677       residual->SetBinContent(i, 0);
678       residual->SetBinError(i, 0);
679     }
680   }
681
682   residual->GetYaxis()->SetTitle("Residuals   1/e (M - R #otimes U)");
683   residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
684   residual->DrawCopy();
685
686   TLine* line = new TLine(-0.5, 0, 150.5, 0);
687   line->SetLineWidth(2);
688   line->Draw();
689
690   pad3->cd();
691   residualHist->SetStats(kFALSE);
692   residualHist->GetXaxis()->SetLabelSize(0.08);
693   residualHist->Fit("gaus");
694   residualHist->Draw();
695
696   canvas->Modified();
697   canvas->SaveAs(canvas->GetName());
698
699   //const char* epsName2 = "proj.eps";
700   //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
701   //canvas->SetGridx();
702   //canvas->SetGridy();
703
704   //canvas->SaveAs(canvas->GetName());
705 }
706
707 void bayesianExample()
708 {
709   TStopwatch watch;
710   watch.Start();
711
712   gSystem->Load("libPWG0base");
713
714   TFile::Open(correctionFile);
715   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
716   mult->LoadHistograms("Multiplicity");
717
718   TFile::Open(measuredFile);
719   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
720   mult2->LoadHistograms("Multiplicity");
721
722   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
723
724   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
725
726   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
727   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
728
729   //mult->DrawComparison("bayesianExample", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
730   DrawResultRatio(mcHist, result, "bayesianExample.eps");
731
732   //Printf("KolmogorovTest says PROB = %f", mcHist->KolmogorovTest(result, "D"));
733   //Printf("Chi2Test says PROB = %f", mcHist->Chi2Test(result));
734
735   // draw residual plot
736
737   // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
738   TH2* convoluted = mult->CalculateMultiplicityESD(result, etaRange);
739   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
740
741   TH1* measured = mult2->GetMultiplicityESD(etaRange)->ProjectionY("measured");
742
743   DrawResiduals(measured, convolutedProj, "bayesianResiduals.eps");
744
745   watch.Stop();
746   watch.Print();
747 }
748
749 void chi2FluctuationResult()
750 {
751   gSystem->Load("libPWG0base");
752
753   TFile::Open(correctionFile);
754   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
755   mult->LoadHistograms("Multiplicity");
756
757   TFile::Open(measuredFile);
758   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
759   mult2->LoadHistograms("Multiplicity");
760
761   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
762   mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
763   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
764
765   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
766   //TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
767
768   mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
769
770   TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_3");
771   canvas->SaveAs("chi2FluctuationResult.eps");
772 }
773
774 void chi2Example()
775 {
776   TStopwatch watch;
777   watch.Start();
778
779   gSystem->Load("libPWG0base");
780
781   TFile::Open(correctionFile);
782   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
783   mult->LoadHistograms("Multiplicity");
784
785   TFile::Open(measuredFile);
786   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
787   mult2->LoadHistograms("Multiplicity");
788
789   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
790   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
791   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
792
793   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
794   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
795
796   DrawResultRatio(mcHist, result, "chi2Example.eps");
797
798   watch.Stop();
799   watch.Print();
800 }
801
802 void chi2ExampleTPC()
803 {
804   TStopwatch watch;
805   watch.Start();
806
807   gSystem->Load("libPWG0base");
808
809   TFile::Open(correctionFileTPC);
810   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
811   mult->LoadHistograms("Multiplicity");
812
813   TFile::Open(measuredFileTPC);
814   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
815   mult2->LoadHistograms("Multiplicity");
816
817   mult->SetMultiplicityESD(etaRangeTPC, mult2->GetMultiplicityESD(etaRangeTPC));
818   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
819   mult->ApplyMinuitFit(etaRangeTPC, kFALSE, AliMultiplicityCorrection::kTrVtx);
820
821   TH1* mcHist = mult2->GetMultiplicityVtx(etaRangeTPC)->ProjectionY("mymc");
822   TH1* result = mult->GetMultiplicityESDCorrected(etaRangeTPC);
823
824   DrawResultRatio(mcHist, result, "chi2ExampleTPC.eps");
825
826   watch.Stop();
827   watch.Print();
828 }
829
830 void bayesianNBD()
831 {
832   gSystem->Load("libPWG0base");
833   TFile::Open("multiplicityMC_3M.root");
834
835   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
836   mult->LoadHistograms("Multiplicity");
837
838   TFile::Open("multiplicityMC_3M_NBD.root");
839   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
840   mult2->LoadHistograms("Multiplicity");
841
842   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
843   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1.0, 100);
844
845   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
846
847   mcHist->Sumw2();
848   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
849
850   //mult->DrawComparison("bayesianNBD", etaRange, kFALSE, kTRUE, mcHist);
851   DrawResultRatio(mcHist, result, "bayesianNBD.eps");
852 }
853
854 void minimizationNBD()
855 {
856   gSystem->Load("libPWG0base");
857   TFile::Open("multiplicityMC_3M.root");
858
859   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
860   mult->LoadHistograms("Multiplicity");
861
862   TFile::Open("multiplicityMC_3M_NBD.root");
863   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
864   mult2->LoadHistograms("Multiplicity");
865
866   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
867   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
868   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
869
870   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
871
872   mcHist->Sumw2();
873   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
874
875   //mult->DrawComparison("minimizationNBD", etaRange, kFALSE, kTRUE, mcHist);
876   DrawResultRatio(mcHist, result, "minimizationNBD.eps");
877 }
878
879 void minimizationInfluenceAlpha()
880 {
881   gSystem->Load("libPWG0base");
882
883   TFile::Open(measuredFile);
884   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
885   mult2->LoadHistograms("Multiplicity");
886
887   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
888   mcHist->Scale(1.0 / mcHist->Integral());
889   mcHist->GetXaxis()->SetRangeUser(0, 200);
890   mcHist->SetStats(kFALSE);
891   mcHist->SetTitle(";true multiplicity;P_{N}");
892
893   TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 300);
894   canvas->Divide(3, 1);
895
896   TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
897
898   TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_100.000000");
899   TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_03_2_100000.000000");
900   TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_06_2_100000000.000000");
901
902   mcHist->Rebin(2);  mcHist->Scale(0.5);
903   hist1->Rebin(2);   hist1->Scale(0.5);
904   hist2->Rebin(2);   hist2->Scale(0.5);
905   hist3->Rebin(2);   hist3->Scale(0.5);
906
907   mcHist->GetXaxis()->SetRangeUser(0, 200);
908
909   canvas->cd(1);
910   gPad->SetLogy();
911   mcHist->Draw();
912   hist1->SetMarkerStyle(5);
913   hist1->SetMarkerColor(2);
914   hist1->Draw("SAME PE");
915
916   canvas->cd(2);
917   gPad->SetLogy();
918   mcHist->Draw();
919   hist2->SetMarkerStyle(5);
920   hist2->SetMarkerColor(2);
921   hist2->Draw("SAME PE");
922
923   canvas->cd(3);
924   gPad->SetLogy();
925   mcHist->Draw();
926   hist3->SetMarkerStyle(5);
927   hist3->SetMarkerColor(2);
928   hist3->Draw("SAME PE");
929
930   canvas->SaveAs("minimizationInfluenceAlpha.eps");
931 }
932
933 void NBDFit()
934 {
935   gSystem->Load("libPWG0base");
936
937   TFile::Open(correctionFile);
938   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
939   mult->LoadHistograms("Multiplicity");
940
941   TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
942   fCurrentESD->Sumw2();
943   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
944
945   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])");
946   func->SetParNames("scaling", "averagen", "k");
947   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
948   func->SetParLimits(1, 0.001, 1000);
949   func->SetParLimits(2, 0.001, 1000);
950   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
951
952   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
953   lognormal->SetParNames("scaling", "mean", "sigma");
954   lognormal->SetParameters(1, 1, 1);
955   lognormal->SetParLimits(0, 0, 10);
956   lognormal->SetParLimits(1, 0, 100);
957   lognormal->SetParLimits(2, 1e-3, 10);
958
959   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
960   fCurrentESD->SetStats(kFALSE);
961   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
962   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
963   fCurrentESD->Draw("HIST");
964   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
965   fCurrentESD->Fit(func, "W0", "", 0, 50);
966   func->SetRange(0, 100);
967   func->Draw("SAME");
968   printf("chi2 = %f\n", func->GetChisquare());
969
970   fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
971   lognormal->SetLineColor(2);
972   lognormal->SetLineStyle(2);
973   lognormal->SetRange(0, 100);
974   lognormal->Draw("SAME");
975
976   canvas->SaveAs("NBDFit.eps");
977 }
978
979 void DifferentSamples()
980 {
981   // data generated by runMultiplicitySelector.C DifferentSamples
982
983   const char* name = "DifferentSamples";
984
985   TFile* file = TFile::Open(Form("%s.root", name));
986
987   TCanvas* canvas = new TCanvas(name, name, 800, 600);
988   canvas->Divide(2, 2);
989
990   for (Int_t i=0; i<4; ++i)
991   {
992     canvas->cd(i+1);
993     gPad->SetTopMargin(0.05);
994     gPad->SetRightMargin(0.05);
995     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", i));
996     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", i));
997     TH1* mc = (TH1*) file->Get(Form("mc_%d", i));
998     mc->Sumw2();
999
1000     chi2Result->Divide(chi2Result, mc, 1, 1, "");
1001     bayesResult->Divide(bayesResult, mc, 1, 1, "");
1002
1003     chi2Result->SetTitle(";true multiplicity;unfolded measured/MC");
1004     chi2Result->GetXaxis()->SetRangeUser(0, 150);
1005     chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1006     chi2Result->GetYaxis()->SetTitleOffset(1.2);
1007     chi2Result->SetLineColor(1);
1008     chi2Result->SetStats(kFALSE);
1009
1010     bayesResult->SetStats(kFALSE);
1011     bayesResult->SetLineColor(2);
1012
1013     chi2Result->DrawCopy("HIST");
1014     bayesResult->DrawCopy("SAME HIST");
1015
1016     TLine* line = new TLine(0, 1, 150, 1);
1017     line->Draw();
1018   }
1019
1020   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1021 }
1022
1023 void StartingConditions()
1024 {
1025   // data generated by runMultiplicitySelector.C StartingConditions
1026
1027   const char* name = "StartingConditions";
1028
1029   TFile* file = TFile::Open(Form("%s.root", name));
1030
1031   TCanvas* canvas = new TCanvas(name, name, 800, 400);
1032   canvas->Divide(2, 1);
1033
1034   TH1* mc = (TH1*) file->Get("mc");
1035   mc->Sumw2();
1036   mc->Scale(1.0 / mc->Integral());
1037
1038   //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
1039
1040   for (Int_t i=0; i<6; ++i)
1041   {
1042     Int_t id = i;
1043     if (id > 2)
1044       id += 2;
1045
1046     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1047     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
1048
1049     chi2Result->Divide(chi2Result, mc, 1, 1, "");
1050     bayesResult->Divide(bayesResult, mc, 1, 1, "");
1051
1052     chi2Result->SetTitle("a) #chi^{2} minimization;true multiplicity;unfolded / MC");
1053     chi2Result->GetXaxis()->SetRangeUser(0, 150);
1054     chi2Result->GetYaxis()->SetRangeUser(0.8, 1.2);
1055     chi2Result->GetYaxis()->SetTitleOffset(1.5);
1056     //chi2Result->SetMarkerStyle(marker[i]);
1057     chi2Result->SetLineColor(i+1);
1058     chi2Result->SetStats(kFALSE);
1059
1060     bayesResult->SetTitle("b) Bayesian method;true multiplicity;unfolded / MC");
1061     bayesResult->GetXaxis()->SetRangeUser(0, 150);
1062     bayesResult->GetYaxis()->SetRangeUser(0.8, 1.2);
1063     bayesResult->GetYaxis()->SetTitleOffset(1.5);
1064     bayesResult->SetStats(kFALSE);
1065     bayesResult->SetLineColor(2);
1066     bayesResult->SetLineColor(i+1);
1067
1068     canvas->cd(1);
1069     gPad->SetLeftMargin(0.12);
1070     chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1071
1072     canvas->cd(2);
1073     gPad->SetLeftMargin(0.12);
1074     bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1075
1076     //TLine* line = new TLine(0, 1, 150, 1);
1077     //line->Draw();
1078   }
1079
1080   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1081 }
1082
1083 void StatisticsPlot()
1084 {
1085   const char* name = "StatisticsPlot";
1086
1087   TFile* file = TFile::Open(Form("%s.root", name));
1088
1089   TCanvas* canvas = new TCanvas(name, name, 600, 400);
1090
1091   TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
1092   fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1093   fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1094   fitResultsChi2->Draw("AP");
1095
1096   TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1097   fitResultsChi2->Fit(f);
1098
1099   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1100
1101   TH1* mc[5];
1102   TH1* result[5];
1103
1104   const char* plotname = "chi2Result";
1105
1106   name = "StatisticsPlotRatios";
1107   canvas = new TCanvas(name, name, 600, 400);
1108
1109   for (Int_t i=0; i<5; ++i)
1110   {
1111     mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1112     result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1113
1114     result[i]->SetLineColor(i+1);
1115     result[i]->Draw(((i == 0) ? "" : "SAME"));
1116   }
1117 }
1118
1119 void SystematicLowEfficiency()
1120 {
1121   gSystem->Load("libPWG0base");
1122
1123   TFile::Open(correctionFile);
1124   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1125   mult->LoadHistograms("Multiplicity");
1126
1127   // calculate result with systematic effect
1128   TFile::Open("multiplicityMC_100k_1_loweff.root");
1129   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1130   mult2->LoadHistograms("Multiplicity");
1131
1132   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1133   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1134   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1135
1136   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1137   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1138
1139   DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
1140
1141   // calculate normal result
1142   TFile::Open("multiplicityMC_100k_1.root");
1143   mult2->LoadHistograms("Multiplicity");
1144
1145   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1146   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1147
1148   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1149
1150   DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
1151
1152   Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
1153 }
1154
1155 void SystematicMisalignment()
1156 {
1157   gSystem->Load("libPWG0base");
1158
1159   TFile::Open(correctionFile);
1160   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1161   mult->LoadHistograms("Multiplicity");
1162
1163   TFile::Open("multiplicityMC_100k_fullmis.root");
1164   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1165   mult2->LoadHistograms("Multiplicity");
1166
1167   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1168   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1169   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1170
1171   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1172   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1173
1174   DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1175 }
1176
1177 void SystematicMisalignmentTPC()
1178 {
1179   gSystem->Load("libPWG0base");
1180
1181   SetTPC();
1182
1183   TFile::Open(correctionFile);
1184   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1185   mult->LoadHistograms("Multiplicity");
1186
1187   TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1188   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1189   mult2->LoadHistograms("Multiplicity");
1190
1191   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1192   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1193   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1194
1195   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1196   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1197
1198   DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1199 }
1200
1201 void EfficiencySpecies()
1202 {
1203   gSystem->Load("libPWG0base");
1204
1205   Int_t marker[] = {24, 25, 26};
1206   Int_t color[] = {1, 2, 4};
1207
1208   // SPD TPC
1209   const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
1210   Float_t etaRange[] = {0.49, 0.9};
1211   const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1212
1213   TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1214   canvas->Divide(2, 1);
1215
1216   for (Int_t loop=0; loop<2; ++loop)
1217   {
1218     Printf("%s", fileName[loop]);
1219
1220     AliCorrection* correction[4];
1221
1222     canvas->cd(loop+1);
1223
1224     gPad->SetGridx();
1225     gPad->SetGridy();
1226     gPad->SetRightMargin(0.05);
1227     //gPad->SetTopMargin(0.05);
1228
1229     TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
1230     legend->SetFillColor(0);
1231     legend->SetEntrySeparation(0.2);
1232
1233     Float_t below = 0;
1234     Float_t total = 0;
1235
1236     TFile* file = TFile::Open(fileName[loop]);
1237     if (!file)
1238     {
1239       Printf("Could not open %s", fileName[loop]);
1240       return;
1241     }
1242
1243     for (Int_t i=0; i<3; ++i)
1244     {
1245       Printf("correction %d", i);
1246
1247       TString name; name.Form("correction_%d", i);
1248       correction[i] = new AliCorrection(name, name);
1249       correction[i]->LoadHistograms();
1250
1251       TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1252       TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1253
1254       // limit vtx axis
1255       gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
1256       meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
1257
1258       // 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)
1259       /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1260         for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1261         {
1262           gene->SetBinContent(x, 0, z, 0);
1263           gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1264           meas->SetBinContent(x, 0, z, 0);
1265           meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1266         }*/
1267
1268       // limit eta axis
1269       gene->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1270       meas->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1271
1272       TH1* genePt = gene->Project3D(Form("z_%d", i));
1273       TH1* measPt = meas->Project3D(Form("z_%d", i));
1274
1275       genePt->Sumw2();
1276       measPt->Sumw2();
1277
1278       TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1279       effPt->Reset();
1280       effPt->Divide(measPt, genePt, 1, 1, "B");
1281
1282       Int_t bin = 0;
1283       for (bin=20; bin>=1; bin--)
1284       {
1285         if (effPt->GetBinContent(bin) < 0.5)
1286           break;
1287       }
1288
1289       Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1290
1291       Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1292       Printf("%.4f of the particles are below that momentum", fraction);
1293
1294       below += genePt->Integral(1, bin);
1295       total += genePt->Integral();
1296
1297       effPt->SetLineColor(color[i]);
1298       effPt->SetMarkerColor(color[i]);
1299       effPt->SetMarkerStyle(marker[i]);
1300
1301       effPt->GetXaxis()->SetRangeUser(0.06, 1);
1302       effPt->GetYaxis()->SetRangeUser(0, 1);
1303
1304       effPt->GetYaxis()->SetTitleOffset(1.2);
1305
1306       effPt->SetStats(kFALSE);
1307       effPt->SetTitle(titles[loop]);
1308       effPt->GetYaxis()->SetTitle("Efficiency");
1309
1310       effPt->DrawCopy((i == 0) ? "" : "SAME");
1311
1312       legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")));
1313     }
1314
1315     Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1316
1317     legend->Draw();
1318   }
1319
1320   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1321 }
1322
1323 void ParticleSpeciesComparison1(Bool_t chi2 = kTRUE, const char* fileNameMC = "multiplicityMC_400k_syst_species.root", const char* fileNameESD = "multiplicityMC_100k_syst.root")
1324 {
1325   gSystem->Load("libPWG0base");
1326
1327   TFile::Open(fileNameESD);
1328   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
1329   TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
1330
1331   TH1* results[10];
1332
1333   // loop over cases (normal, enhanced/reduced ratios)
1334   Int_t nMax = 7;
1335   for (Int_t i = 0; i<nMax; ++i)
1336   {
1337     TString folder;
1338     folder.Form("Multiplicity_%d", i);
1339
1340     AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
1341
1342     TFile::Open(fileNameMC);
1343     mult->LoadHistograms();
1344
1345     mult->SetMultiplicityESD(etaRange, hist);
1346
1347     if (chi2)
1348     {
1349       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1350       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1351       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1352     }
1353     else
1354     {
1355       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1356       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1357     }
1358
1359     //Float_t averageRatio = 0;
1360     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1361
1362     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1363
1364     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1365   }
1366
1367   DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
1368
1369   TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
1370
1371   for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1372   {
1373     results[0]->SetBinError(i, 0);
1374     mc->SetBinError(i, 0);
1375   }
1376
1377   const char* legendStrings[] = { "#pi^{#pm}", 0, "K^{#pm}", 0, "p,#bar{p}", 0 };
1378
1379   DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1380
1381   //not valid: draw chi2 uncertainty on top!
1382   /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1383   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1384   errorHist->SetLineColor(1);
1385   errorHist->SetLineWidth(2);
1386   TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1387   for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1388   {
1389     errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1390     errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1391   }
1392
1393   errorHist->DrawCopy("SAME");
1394   errorHist2->DrawCopy("SAME");*/
1395
1396   //canvas->SaveAs(canvas->GetName());
1397
1398   DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1399
1400   //errorHist->DrawCopy("SAME");
1401   //errorHist2->DrawCopy("SAME");
1402
1403   //canvas2->SaveAs(canvas2->GetName());
1404 }
1405
1406 /*void ParticleSpeciesComparison2()
1407 {
1408   gSystem->Load("libPWG0base");
1409
1410   const char* fileNameMC = "multiplicityMC_400k_syst.root";
1411   const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1412   Bool_t chi2 = 0;
1413
1414   TFile::Open(fileNameMC);
1415   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1416   mult->LoadHistograms();
1417
1418   TH1* mc[10];
1419   TH1* results[10];
1420
1421   // loop over cases (normal, enhanced/reduced ratios)
1422   Int_t nMax = 7;
1423   for (Int_t i = 0; i<nMax; ++i)
1424   {
1425     TString folder;
1426     folder.Form("Multiplicity_%d", i);
1427
1428     AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1429
1430     TFile::Open(fileNameESD);
1431     mult2->LoadHistograms();
1432
1433     mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1434
1435     if (chi2)
1436     {
1437       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1438       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1439       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1440     }
1441     else
1442     {
1443       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1444       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1445     }
1446
1447     //Float_t averageRatio = 0;
1448     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1449
1450     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1451
1452     TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1453     mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1454
1455     //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1456     //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1457
1458     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1459   }
1460
1461   DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1462 }*/
1463
1464 TH1* Invert(TH1* eff)
1465 {
1466   // calculate corr = 1 / eff
1467
1468   TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1469   corr->Reset();
1470
1471   for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1472   {
1473     if (eff->GetBinContent(i) > 0)
1474     {
1475       corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1476       corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1477     }
1478   }
1479
1480   return corr;
1481 }
1482
1483 void TriggerVertexCorrection()
1484 {
1485   //
1486   // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1487   //
1488
1489   gSystem->Load("libPWG0base");
1490
1491   TFile::Open(correctionFile);
1492   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1493   mult->LoadHistograms("Multiplicity");
1494
1495   TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1496   TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1497
1498   TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
1499
1500   corrINEL->SetStats(kFALSE);
1501   corrINEL->GetXaxis()->SetRangeUser(0, 30);
1502   corrINEL->GetYaxis()->SetRangeUser(0, 5);
1503   corrINEL->SetTitle(";true multiplicity;correction factor");
1504   corrINEL->SetMarkerStyle(22);
1505   corrINEL->Draw("PE");
1506
1507   corrMB->SetStats(kFALSE);
1508   corrMB->SetLineColor(2);
1509   corrMB->SetMarkerStyle(25);
1510   corrMB->SetMarkerColor(2);
1511   corrMB->Draw("SAME PE");
1512
1513   TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
1514   legend->SetFillColor(0);
1515   legend->AddEntry(corrINEL, "correction to inelastic sample");
1516   legend->AddEntry(corrMB, "correction to minimum bias sample");
1517
1518   legend->Draw();
1519
1520   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1521 }
1522
1523 void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1524 {
1525   gSystem->Load("libPWG0base");
1526
1527   TFile::Open(correctionFile);
1528   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1529   mult->LoadHistograms("Multiplicity");
1530
1531   TFile::Open(measuredFile);
1532   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1533   mult2->LoadHistograms("Multiplicity");
1534
1535   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1536
1537   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1538
1539   TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1540
1541   TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1542   TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1543
1544   if (!mc)
1545   {
1546     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1547     DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1548   }
1549
1550   TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1551   canvas->SetGridx();
1552   canvas->SetGridy();
1553   canvas->SetRightMargin(0.05);
1554   canvas->SetTopMargin(0.05);
1555
1556   errorResponse->SetLineColor(1);
1557   errorResponse->GetXaxis()->SetRangeUser(0, 200);
1558   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1559   errorResponse->SetStats(kFALSE);
1560   errorResponse->SetTitle(";true multiplicity;Uncertainty");
1561
1562   errorResponse->Draw();
1563
1564   errorMeasured->SetLineColor(2);
1565   errorMeasured->Draw("SAME");
1566
1567   errorBoth->SetLineColor(4);
1568   errorBoth->Draw("SAME");
1569
1570   Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
1571   Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
1572   Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
1573
1574   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1575
1576   TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
1577   errorResponse->Write();
1578   errorMeasured->Write();
1579   errorBoth->Write();
1580   file->Close();
1581 }
1582
1583 void StatisticalUncertaintyCompare(const char* det = "SPD")
1584 {
1585   TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1586   TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1587   TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1588   TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1589
1590   TString str;
1591   str.Form("StatisticalUncertaintyCompare%s", det);
1592
1593   TCanvas* canvas = new TCanvas(str, str, 600, 400);
1594   canvas->SetGridx();
1595   canvas->SetGridy();
1596   canvas->SetRightMargin(0.05);
1597   canvas->SetTopMargin(0.05);
1598
1599   errorResponse->SetLineColor(1);
1600   errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
1601   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1602   errorResponse->SetStats(kFALSE);
1603   errorResponse->GetYaxis()->SetTitleOffset(1.2);
1604   errorResponse->SetTitle(";true multiplicity;#sigma(U-T)/T");
1605
1606   errorResponse->Draw();
1607
1608   errorMeasured->SetLineColor(2);
1609   errorMeasured->Draw("SAME");
1610
1611   errorBoth->SetLineColor(4);
1612   errorBoth->Draw("SAME");
1613
1614   TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
1615   TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
1616
1617   errorBoth2->SetLineColor(4);
1618   errorBoth2->SetLineStyle(2);
1619   errorBoth2->Draw("SAME");
1620
1621   TLegend* legend = new TLegend(0.2, 0.6, 0.6, 0.9);
1622   legend->SetFillColor(0);
1623   legend->AddEntry(errorResponse, "response matrix (Bayesian)");
1624   legend->AddEntry(errorMeasured, "measured (Bayesian)");
1625   legend->AddEntry(errorBoth, "both (Bayesian)");
1626   legend->AddEntry(errorBoth2, "both (#chi^{2} minimization)");
1627   legend->Draw();
1628
1629   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1630 }
1631
1632 void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
1633 {
1634  const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
1635
1636   gSystem->Load("libPWG0base");
1637
1638   TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
1639   canvas->SetGridx();
1640   canvas->SetGridy();
1641   canvas->SetRightMargin(0.05);
1642   canvas->SetTopMargin(0.05);
1643
1644   AliMultiplicityCorrection* data[4];
1645   TH1* effArray[4];
1646
1647   Int_t markers[] = { 24, 25, 26, 5 };
1648   Int_t colors[] = { 1, 2, 3, 4 };
1649
1650   TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
1651   legend->SetFillColor(0);
1652
1653   TH1* effError = 0;
1654
1655   for (Int_t i=0; i<4; ++i)
1656   {
1657     TString name;
1658     name.Form("Multiplicity_%d", i);
1659
1660     TFile::Open(files[i]);
1661     data[i] = new AliMultiplicityCorrection(name, name);
1662
1663     if (i < 3)
1664     {
1665       data[i]->LoadHistograms("Multiplicity");
1666     }
1667     else
1668       data[i]->LoadHistograms("Multiplicity_0");
1669
1670     TH1* eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
1671     effArray[i] = eff;
1672
1673     eff->GetXaxis()->SetRangeUser(0, 15);
1674     eff->GetYaxis()->SetRangeUser(0, 1.1);
1675     eff->SetStats(kFALSE);
1676     eff->SetTitle(";true multiplicity;Efficiency");
1677     eff->SetLineColor(colors[i]);
1678     eff->SetMarkerColor(colors[i]);
1679     eff->SetMarkerStyle(markers[i]);
1680
1681     if (i == 3)
1682     {
1683       for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1684         eff->SetBinError(bin, 0);
1685
1686       // loop over cross section combinations
1687       for (Int_t j=1; j<7; ++j)
1688       {
1689         AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
1690         mult->LoadHistograms(Form("Multiplicity_%d", j));
1691
1692         TH1* eff2 = mult->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType);
1693
1694         for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1695         {
1696           // TODO we could also do asymmetric errors here
1697           Float_t deviation = TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
1698
1699           eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), (Double_t) deviation));
1700         }
1701       }
1702
1703       for (Int_t bin=1; bin<=20; bin++)
1704         if (eff->GetBinContent(bin) > 0)
1705           Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1706       
1707       if (uncertainty) {
1708         effError = (TH1*) eff->Clone("effError");
1709         effError->Reset();
1710
1711         for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
1712           if (eff->GetBinContent(bin) > 0)
1713             effError->SetBinContent(bin, 10.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1714
1715         effError->SetLineColor(1);
1716         effError->SetMarkerStyle(1);
1717         effError->DrawCopy("SAME HIST");
1718       }
1719     }
1720
1721     eff->SetBinContent(1, 0);
1722     eff->SetBinError(1, 0);
1723
1724     canvas->cd();
1725     if (i == 0)
1726     {
1727       eff->DrawCopy("P");
1728     }
1729     else
1730       eff->DrawCopy("SAME P");
1731
1732     legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
1733   }
1734
1735   if (uncertainty)
1736     legend->AddEntry(effError, "relative syst. uncertainty #times 10");
1737
1738   legend->Draw();
1739
1740   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1741 }
1742
1743 void ModelDependencyPlot()
1744 {
1745   gSystem->Load("libPWG0base");
1746
1747   TFile::Open("multiplicityMC_3M.root");
1748   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1749   mult->LoadHistograms("Multiplicity");
1750
1751   TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
1752
1753   TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
1754   canvas->SetGridx();
1755   canvas->SetGridy();
1756   //canvas->SetRightMargin(0.05);
1757   //canvas->SetTopMargin(0.05);
1758
1759   canvas->Divide(2, 1);
1760
1761   canvas->cd(2);
1762   gPad->SetLogy();
1763  
1764   Int_t selectedMult = 30;
1765   Int_t yMax = 200000;
1766
1767   TH1* full = proj->ProjectionX("full");
1768   TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult)); 
1769
1770   full->SetStats(kFALSE);
1771   full->GetXaxis()->SetRangeUser(0, 200);
1772   full->GetYaxis()->SetRangeUser(5, yMax);
1773   full->SetTitle(";multiplicity");
1774
1775   selected->SetLineColor(0);
1776   selected->SetMarkerColor(2);
1777   selected->SetMarkerStyle(7);
1778
1779   full->Draw();
1780   selected->Draw("SAME P");
1781
1782   TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
1783   legend->SetFillColor(0);
1784   legend->AddEntry(full, "true");
1785   legend->AddEntry(selected, "measured");
1786   legend->Draw();
1787  
1788   TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
1789   line->SetLineWidth(2);
1790   line->Draw();
1791
1792   canvas->cd(1);
1793   gPad->SetLogy();
1794
1795   full = proj->ProjectionY("full2");
1796   selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
1797
1798   full->SetStats(kFALSE);
1799   full->GetXaxis()->SetRangeUser(0, 200);
1800   full->GetYaxis()->SetRangeUser(5, yMax);
1801   full->SetTitle(";multiplicity");
1802
1803   full->SetLineColor(0);
1804   full->SetMarkerColor(2);
1805   full->SetMarkerStyle(7);
1806
1807   full->Draw("P");
1808   selected->Draw("SAME");
1809
1810   legend->Draw();
1811
1812   line = new TLine(selectedMult, 5, selectedMult, yMax);
1813   line->SetLineWidth(2);
1814   line->Draw();
1815
1816   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1817 }
1818
1819 void SystematicpTSpectrum()
1820 {
1821   gSystem->Load("libPWG0base");
1822
1823   TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
1824   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1825   mult->LoadHistograms("Multiplicity");
1826
1827   TFile::Open("multiplicityMC_100k_syst.root");
1828   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1829   mult2->LoadHistograms("Multiplicity");
1830
1831   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1832   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1833   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1834
1835   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1836   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1837
1838   DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
1839 }
1840
1841 // to be deleted
1842 /*void covMatrix(Bool_t mc = kTRUE)
1843 {
1844   gSystem->Load("libPWG0base");
1845
1846   TFile::Open(correctionFile);
1847   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1848   mult->LoadHistograms("Multiplicity");
1849
1850   TFile::Open(measuredFile);
1851   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1852   mult2->LoadHistograms("Multiplicity");
1853
1854   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1855
1856   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1857
1858   mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
1859 }*/
1860
1861 Double_t FitPtFunc(Double_t *x, Double_t *par)
1862 {
1863   Double_t xx = x[0];
1864
1865   Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
1866   Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
1867
1868   const Float_t kTransitionWidth = 0;
1869
1870   // power law part
1871   if (xx < par[0] - kTransitionWidth)
1872   {
1873     return val1;
1874   }
1875   /*else if (xx < par[0] + kTransitionWidth)
1876   {
1877     // smooth transition
1878     Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
1879     return (1 - factor) * val1 + factor * val2;
1880   }*/
1881   else
1882   {
1883     return val2;
1884   }
1885 }
1886
1887 void FitPt(const char* fileName = "firstplots100k_truept.root")
1888 {
1889   gSystem->Load("libPWG0base");
1890
1891   TFile::Open(fileName);
1892
1893   /*
1894   // merge corrections
1895   AliCorrection* correction[4];
1896   TList list;
1897
1898   for (Int_t i=0; i<4; ++i)
1899   {
1900     Printf("correction %d", i);
1901
1902     TString name; name.Form("correction_%d", i);
1903     correction[i] = new AliCorrection(name, name);
1904     correction[i]->LoadHistograms();
1905
1906     if (i > 0)
1907       list.Add(correction[i]);
1908   }
1909
1910   correction[0]->Merge(&list);
1911
1912   TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
1913
1914   // limit vtx, eta axis
1915   gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
1916   gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
1917
1918   TH1* genePt = gene->Project3D("z");*/
1919   TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
1920   genePt->Sumw2();
1921
1922   //genePt->Scale(1.0 / genePt->Integral());
1923
1924   // normalize by bin width
1925   for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1926     genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1927
1928   /// genePt->GetXaxis()->GetBinCenter(x));
1929
1930   genePt->GetXaxis()->SetRangeUser(0, 7.9);
1931   //genePt->GetYaxis()->SetTitle("a.u.");
1932
1933   TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
1934   //func->SetLineColor(2);
1935   func->SetParameters(1, -1, 1, 1, 1);
1936   func->SetParLimits(3, 1, 10);
1937   func->SetParLimits(4, 0, 10);
1938
1939   //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
1940
1941   //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
1942   //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
1943   //func->FixParameter(0, 0.314);
1944   //func->SetParLimits(0, 0.1, 0.3);
1945
1946   genePt->SetMarkerStyle(25);
1947   genePt->SetTitle("");
1948   genePt->SetStats(kFALSE);
1949   genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
1950   //func->Draw("SAME");
1951
1952   // fit only exp. part
1953   func->SetParameters(1, -1);
1954   func->FixParameter(2, 0);
1955   func->FixParameter(3, 1);
1956   func->FixParameter(4, 1);
1957   genePt->Fit(func, "0", "", 0.2, 1);
1958
1959   new TCanvas;
1960   genePt->DrawCopy("P");
1961   func->SetRange(0.02, 8);
1962   func->DrawCopy("SAME");
1963   gPad->SetLogy();
1964
1965   // now fix exp. parameters and fit second part
1966   Double_t param0 = func->GetParameter(0);
1967   Double_t param1 = func->GetParameter(1);
1968   func->SetParameters(0, -1, 1, 1, 1);
1969   func->FixParameter(0, 0);
1970   func->FixParameter(1, -1);
1971   func->ReleaseParameter(2);
1972   func->ReleaseParameter(3);
1973   func->ReleaseParameter(4);
1974   func->SetParLimits(3, 1, 10);
1975   func->SetParLimits(4, 0, 10);
1976
1977   genePt->Fit(func, "0", "", 1.5, 4);
1978
1979   new TCanvas;
1980   genePt->DrawCopy("P");
1981   func->SetRange(0.02, 8);
1982   func->DrawCopy("SAME");
1983   gPad->SetLogy();
1984
1985   // fit both
1986   func->SetParameter(0, param0);
1987   func->SetParameter(1, param1);
1988   func->ReleaseParameter(0);
1989   func->ReleaseParameter(1);
1990
1991   new TCanvas;
1992   genePt->DrawCopy("P");
1993   func->SetRange(0.02, 5);
1994   func->DrawCopy("SAME");
1995   gPad->SetLogy();
1996
1997   genePt->Fit(func, "0", "", 0.2, 4);
1998
1999   TCanvas* canvas = new TCanvas("FitPt", "FitPt", 800, 400);
2000   canvas->Divide(2, 1);
2001   canvas->cd(1);
2002
2003   gPad->SetGridx();
2004   gPad->SetGridy();
2005   gPad->SetLeftMargin(0.13);
2006   gPad->SetRightMargin(0.05);
2007   gPad->SetTopMargin(0.05);
2008
2009   genePt->GetXaxis()->SetRangeUser(0, 4.9);
2010   genePt->GetYaxis()->SetRangeUser(1e-2, 1e4);
2011   genePt->GetYaxis()->SetTitleOffset(1.4);
2012   genePt->GetXaxis()->SetTitleOffset(1.1);
2013   genePt->DrawCopy("P");
2014   func->SetRange(0.02, 5);
2015   func->DrawCopy("SAME");
2016   gPad->SetLogy();
2017
2018   canvas->cd(2);
2019
2020   TH1* genePtClone = (TH1*) genePt->Clone("genePtClone");
2021   genePtClone->Reset();
2022   genePtClone->DrawCopy("P");
2023
2024   gPad->SetGridx();
2025   gPad->SetGridy();
2026   gPad->SetLeftMargin(0.13);
2027   gPad->SetRightMargin(0.05);
2028   gPad->SetTopMargin(0.05);
2029
2030   func->DrawCopy("SAME");
2031   gPad->SetLogy();
2032
2033   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2034
2035   TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2036
2037   TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2038
2039   TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2040
2041   for (Int_t param=0; param<5; param++)
2042   {
2043     for (Int_t sign=0; sign<2; sign++)
2044     {
2045       TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign));  // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2046       func2->SetParameters(func->GetParameters());
2047       //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2048
2049       Float_t factor = ((sign == 0) ? 0.9 : 1.1);
2050       func2->SetParameter(param, func2->GetParameter(param) * factor);
2051       //func2->Print();
2052
2053       canvas->cd(2);
2054       func2->SetLineWidth(1);
2055       func2->SetLineColor(2);
2056       func2->DrawCopy("SAME");
2057
2058       canvas2->cd();
2059       TH1* second = func2->GetHistogram();
2060       second->Divide(first);
2061       second->SetLineColor(param + 1);
2062       second->GetYaxis()->SetRangeUser(0, 2);
2063       second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2064       second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2065     }
2066   }
2067
2068   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2069   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2070
2071   file->Close();
2072 }
2073
2074 void DrawSystematicpT()
2075 {
2076   TFile* file = TFile::Open("SystematicpT.root");
2077
2078   TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2079   TH1* result2 = (TH1*) file->Get("result_unity");
2080
2081   TH1* mcHist[12];
2082   TH1* result[12];
2083
2084   Int_t nParams = 5;
2085
2086   for (Int_t id=0; id<nParams*2; ++id)
2087   {
2088     mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2089     result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2090   }
2091
2092   DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2093
2094   //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2095
2096   DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2097
2098   //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2099
2100   // does not make sense: mc is different
2101   //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2102 }
2103
2104 void SystematicpT(Bool_t chi2 = 1)
2105 {
2106   gSystem->Load("libPWG0base");
2107
2108   TFile::Open("ptspectrum900.root");
2109   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2110   mult->LoadHistograms("Multiplicity");
2111
2112   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2113
2114   TH1* mcHist[12];
2115   TH1* result[12];
2116
2117   Int_t nParams = 5;
2118
2119   for (Int_t param=0; param<nParams; param++)
2120   {
2121     for (Int_t sign=0; sign<2; sign++)
2122     {
2123       // calculate result with systematic effect
2124       TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2125       mult2->LoadHistograms("Multiplicity");
2126
2127       mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2128
2129       if (chi2)
2130       {
2131         mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2132         mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2133       }
2134       else
2135         mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2136
2137       Int_t id = param * 2 + sign;
2138
2139       mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2140       result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2141
2142       TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2143       DrawResultRatio(mcHist[id], result[id], tmp);
2144     }
2145   }
2146
2147   // calculate normal result
2148   TFile::Open("ptspectrum100_1.root");
2149   mult2->LoadHistograms("Multiplicity");
2150   TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2151
2152   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2153
2154   if (chi2)
2155   {
2156     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2157   }
2158   else
2159     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2160
2161   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2162
2163   TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2164   mcHist2->Write();
2165   result2->Write();
2166   for (Int_t id=0; id<nParams*2; ++id)
2167   {
2168     mcHist[id]->Write();
2169     result[id]->Write();
2170   }
2171   file->Close();
2172
2173   DrawSystematicpT();
2174 }
2175
2176 void DrawSystematicpT2()
2177 {
2178   //displayRange = 200;
2179
2180   // read from file
2181   TFile* file = TFile::Open("SystematicpT2.root");
2182   TH1* mcHist = (TH1*) file->Get("mymc");
2183   TH1* result[12];
2184   result[0] = (TH1*) file->Get("result_unity");
2185   Int_t nParams = 5;
2186   for (Int_t id=0; id<nParams*2; ++id)
2187     result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2188
2189   DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2190   DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2191   DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2192 }
2193
2194 void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2195 {
2196   gSystem->Load("libPWG0base");
2197
2198   if (tpc)
2199   {
2200     SetTPC();
2201     TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2202   }
2203   else
2204     TFile::Open("ptspectrum100_1.root");
2205
2206   AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2207   measured->LoadHistograms("Multiplicity");
2208   TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2209
2210   TH1* result[12];
2211
2212   Int_t nParams = 5;
2213
2214   // -1 = unity change, 0...4 parameters
2215   for (Int_t id=-1; id<nParams*2; id++)
2216   {
2217     Int_t param = id / 2;
2218     Int_t sign = id % 2;
2219
2220     TString idStr;
2221     if (id == -1)
2222     {
2223       idStr = "unity";
2224     }
2225     else
2226       idStr.Form("%d_%d", param, sign);
2227
2228     // calculate result with systematic effect
2229     if (tpc)
2230     {
2231       TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2232     }
2233     else
2234       TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2235
2236     AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2237     response->LoadHistograms("Multiplicity");
2238
2239     response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2240
2241     if (chi2)
2242     {
2243       response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2244       response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2245     }
2246     else
2247       response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2248
2249     result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2250
2251     TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2252     DrawResultRatio(mcHist, result[id+1], tmp);
2253   }
2254
2255   TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2256   mcHist->Write();
2257   for (Int_t id=0; id<nParams*2+1; ++id)
2258     result[id]->Write();
2259   file->Close();
2260
2261   DrawSystematicpT2();
2262 }
2263
2264 void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2265 {
2266   // only needed for TPC
2267   SetTPC();
2268
2269   gSystem->Load("libPWG0base");
2270
2271   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2272   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2273   mult->LoadHistograms("Multiplicity");
2274
2275   TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2276   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2277   mult2->LoadHistograms("Multiplicity");
2278
2279   // "normal" result
2280   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2281
2282   if (chi2)
2283   {
2284     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2285     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2286   }
2287   else
2288     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2289
2290   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2291   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2292
2293   TH1* syst[2];
2294
2295   // change of pt spectrum (down)
2296   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2297   AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2298   mult3->LoadHistograms("Multiplicity");
2299   mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2300   if (chi2)
2301   {
2302     mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2303     mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2304   }
2305   else
2306     mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2307   syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2308
2309   // change of pt spectrum (up)
2310   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2311   AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2312   mult4->LoadHistograms("Multiplicity");
2313   mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2314   if (chi2)
2315   {
2316     mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2317     mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2318   }
2319   else
2320     mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2321   syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
2322
2323   DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
2324
2325   Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2326   Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2327 }
2328
2329 TH1* SystematicsSummary(Bool_t tpc = 1)
2330 {
2331   Int_t nEffects = 7;
2332
2333   TH1* effects[10];
2334   const char** names = 0;
2335   Int_t colors[] = { 1, 2, 3, 4, 6, 7, 8 };
2336   Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
2337
2338   for (Int_t i=0; i<nEffects; ++i)
2339     effects[i] = new TH1F("SystematicsSummary", ";true multiplicity;Effect", 201, -0.5, 200.5);
2340
2341   if (tpc)
2342   {
2343     SetTPC();
2344
2345     const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2346     names = namesTPC;
2347
2348     // method
2349     TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2350     TH1* hist = (TH1*) file->Get("errorBoth");
2351
2352     // smooth a bit, but skip 0 bin
2353     effects[0]->SetBinContent(2, hist->GetBinContent(2));
2354     for (Int_t i=3; i<=200; ++i)
2355       effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2356
2357     // relative x-section
2358     effects[1]->SetBinContent(2, 0.005);
2359     effects[1]->SetBinContent(3, 0.0025);
2360     effects[1]->SetBinContent(4, 0.0025);
2361
2362     // particle composition
2363     for (Int_t i=2; i<=101; ++i)
2364     {
2365       if (i < 41)
2366       {
2367         effects[2]->SetBinContent(i, 0.01);
2368       }
2369       else if (i < 76)
2370       {
2371         effects[2]->SetBinContent(i, 0.02);
2372       }
2373       else
2374         effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2375     }
2376
2377     // pt cut off (only tpc)
2378     for (Int_t i=2; i<=101; ++i)
2379     {
2380       if (i < 11)
2381       {
2382         effects[3]->SetBinContent(i, 0.05);
2383       }
2384       else if (i < 51)
2385       {
2386         effects[3]->SetBinContent(i, 0.01);
2387       }
2388       else
2389         effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2390     }
2391
2392     // track selection (only tpc)
2393     for (Int_t i=2; i<=101; ++i)
2394       effects[4]->SetBinContent(i, 0.03);
2395
2396     // secondaries
2397     for (Int_t i=2; i<=101; ++i)
2398       effects[5]->SetBinContent(i, 0.01);
2399
2400     // pt spectrum
2401     for (Int_t i=2; i<=101; ++i)
2402     {
2403       if (i < 21)
2404       {
2405         effects[6]->SetBinContent(i, 0.05);
2406       }
2407       else if (i < 51)
2408       {
2409         effects[6]->SetBinContent(i, 0.02);
2410       }
2411       else
2412         effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2413     }
2414
2415   }
2416   else
2417   {
2418     displayRange = 200;
2419     nEffects = 5;
2420
2421     const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Rel. cross-section", "Particle composition", "Secondaries", "p_{t} spectrum"};
2422     names = namesSPD;
2423
2424     // method
2425     TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2426     TH1* hist = (TH1*) file->Get("errorBoth");
2427
2428     // smooth a bit, but skip 0 bin
2429     effects[0]->SetBinContent(2, hist->GetBinContent(2));
2430     for (Int_t i=3; i<=201; ++i)
2431       effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2432
2433     // relative x-section
2434     effects[1]->SetBinContent(2, 0.01);
2435     effects[1]->SetBinContent(3, 0.005);
2436
2437     // particle composition
2438     for (Int_t i=2; i<=201; ++i)
2439     {
2440       if (i < 6)
2441       {
2442         effects[2]->SetBinContent(i, 0.3);
2443       }
2444       else if (i < 11)
2445       {
2446         effects[2]->SetBinContent(i, 0.05);
2447       }
2448       else if (i < 121)
2449       {
2450         effects[2]->SetBinContent(i, 0.02);
2451       }
2452       else if (i < 151)
2453       {
2454         effects[2]->SetBinContent(i, 0.02 + 0.04 / 30 * (i - 121));
2455       }
2456       else
2457         effects[2]->SetBinContent(i, 0.06 + 0.1 / 30 * (i - 151));
2458     }
2459
2460     // secondaries
2461     for (Int_t i=2; i<=201; ++i)
2462       effects[3]->SetBinContent(i, 0.01);
2463
2464     // pt spectrum
2465     for (Int_t i=2; i<=201; ++i)
2466     {
2467       if (i < 6)
2468       {
2469         effects[4]->SetBinContent(i, 1);
2470       }
2471       else if (i < 121)
2472       {
2473         effects[4]->SetBinContent(i, 0.03);
2474       }
2475       else if (i < 151)
2476       {
2477         effects[4]->SetBinContent(i, 0.03 + 0.07 / 30 * (i - 121));
2478       }
2479       else
2480         effects[4]->SetBinContent(i, 0.1);
2481     }
2482   }
2483
2484   TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 400);
2485   canvas->SetRightMargin(0.25);
2486   canvas->SetTopMargin(0.05);
2487   TLegend* legend = new TLegend(0.2, 0.4, 0.5, 0.4 + 0.5 * nEffects / 7);
2488   legend->SetFillColor(0);
2489
2490   for (Int_t i=0; i<nEffects; ++i)
2491   {
2492     TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2493     /*current->Reset();
2494     for (Int_t j=0; j<nEffects-i; ++j)
2495       current->Add(effects[j]);*/
2496
2497     current->SetLineColor(colors[i]);
2498     //current->SetFillColor(colors[i]);
2499     current->SetMarkerColor(colors[i]);
2500     //current->SetMarkerStyle(markers[i]);
2501
2502     current->SetStats(kFALSE);
2503     current->GetYaxis()->SetRangeUser(0, 0.4);
2504     current->GetXaxis()->SetRangeUser(0, displayRange);
2505     current->DrawCopy(((i == 0) ? "" : "SAME"));
2506     legend->AddEntry(current, names[i]);
2507
2508     TLatex* text = new TLatex(displayRange+5, current->GetBinContent(displayRange+1), names[i]);
2509     text->SetTextColor(colors[i]);
2510     text->Draw();
2511   }
2512
2513   // add total in square
2514   TH1* total = (TH1*) effects[0]->Clone("total");
2515   total->Reset();
2516
2517   for (Int_t i=0; i<nEffects; ++i)
2518   {
2519     //Printf("%d %f", i, effects[i]->GetBinContent(20));
2520     effects[i]->Multiply(effects[i]);
2521     total->Add(effects[i]);
2522   }
2523
2524   for (Int_t i=1; i<=total->GetNbinsX(); ++i)
2525     if (total->GetBinContent(i) > 0)
2526       total->SetBinContent(i, TMath::Min(sqrt(total->GetBinContent(i)), 1.0));
2527
2528   //Printf("%f", total->GetBinContent(20));
2529
2530   total->SetMarkerStyle(3);
2531   total->SetMarkerColor(1);
2532   legend->AddEntry(total, "total");
2533   total->DrawCopy("SAME P");
2534
2535   legend->Draw();
2536
2537   canvas->SaveAs(canvas->GetName());
2538
2539   return total;
2540 }
2541
2542 void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE, Bool_t small = kFALSE)
2543 {
2544   gSystem->Load("libPWG0base");
2545
2546   if (tpc)
2547     SetTPC();
2548
2549   if (!chi2)
2550     Printf("WARNING: Bayesian set. This is only for test!");
2551
2552   // systematic error
2553   TH1* error = SystematicsSummary(tpc);
2554
2555   TFile::Open(correctionFile);
2556   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2557   mult->LoadHistograms("Multiplicity");
2558
2559   TFile::Open(measuredFile);
2560   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2561   mult2->LoadHistograms("Multiplicity");
2562
2563   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2564
2565   if (chi2)
2566   {
2567     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2568     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kINEL);
2569   }
2570   else
2571     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kINEL, 1, 100, 0, kFALSE);
2572
2573   TH1* mcHist = mult2->GetMultiplicityINEL(etaRange)->ProjectionY("mymc");
2574   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2575
2576   DrawResultRatio(mcHist, result, "finalPlotCheck.eps");
2577
2578   // normalize result
2579   result->Scale(1.0 / result->Integral(2, 200));
2580
2581   result->GetXaxis()->SetRangeUser(0, ((tpc) ? displayRange : 200));
2582   result->SetBinContent(1, 0); result->SetBinError(1, 0);
2583   result->SetTitle(";true multiplicity;Probability");
2584   result->SetLineColor(1);
2585   result->SetStats(kFALSE);
2586
2587   TH1* systError = (TH1*) result->Clone("systError");
2588   for (Int_t i=2; i<=systError->GetNbinsX(); ++i)
2589     systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2590
2591   // change error drawing style
2592   systError->SetFillColor(15);
2593
2594   TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", (small) ? 600 : 800, 400);
2595   canvas->SetRightMargin(0.05);
2596   canvas->SetTopMargin(0.05);
2597
2598   systError->Draw("E2 ][");
2599   result->DrawCopy("SAME E ][");
2600   canvas->SetLogy();
2601
2602   //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
2603   TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
2604   text->SetFillColor(0);
2605   text->SetTextAlign(12);
2606   text->AddText("Systematic errors summed quadratically");
2607   text->AddText("0.6 million minimum bias events");
2608   text->AddText("corrected to inelastic events");
2609   text->Draw("B");
2610
2611   TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
2612   text2->SetFillColor(0);
2613   text2->SetTextAlign(12);
2614   text2->AddText("#sqrt{s} = 14 TeV");
2615   if (tpc)
2616   {
2617     text2->AddText("|#eta| < 0.9");
2618   }
2619   else
2620     text2->AddText("|#eta| < 2.0");
2621   text2->AddText("simulated data (PYTHIA)");
2622   text2->Draw("B");
2623
2624   if (tpc)
2625   {
2626     TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
2627     text3->SetNDC();
2628     text3->Draw();
2629   }
2630   else
2631   {
2632     TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
2633     text3->SetNDC();
2634     text3->Draw();
2635   }
2636
2637   // alice logo
2638   TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
2639   pad->Draw();
2640   pad->cd();
2641   TImage* img = TImage::Open("$HOME/alice.png");
2642   img->SetImageQuality(TAttImage::kImgBest);
2643   img->Draw();
2644
2645   canvas->Modified();
2646
2647 /*  TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
2648   text->SetTextSize(0.04);
2649   text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
2650   text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
2651   text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
2652   text->Draw();*/
2653
2654
2655   canvas->SaveAs(canvas->GetName());
2656 }
2657
2658 void BlobelUnfoldingExample()
2659 {
2660   const Int_t kSize = 20;
2661
2662   TMatrixD matrix(kSize, kSize);
2663   for (Int_t x=0; x<kSize; x++)
2664   {
2665     for (Int_t y=0; y<kSize; y++)
2666     {
2667       if (x == y)
2668       {
2669         if (x == 0 || x == kSize -1)
2670         {
2671           matrix(x, y) = 0.75;
2672         }
2673         else
2674           matrix(x, y) = 0.5;
2675       }
2676       else if (TMath::Abs(x - y) == 1)
2677       {
2678         matrix(x, y) = 0.25;
2679       }
2680     }
2681   }
2682
2683   //matrix.Print();
2684
2685   TMatrixD inverted(matrix);
2686   inverted.Invert();
2687
2688   //inverted.Print();
2689
2690   TH1F* inputDist = new TH1F("inputDist", ";t;#tilde{T}(t)", kSize, -0.5, (Float_t) kSize - 0.5);
2691   TVectorD inputDistVector(kSize);
2692   TH1F* unfolded = inputDist->Clone("unfolded");
2693   TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
2694   measuredIdealDist->SetTitle(";m;#tilde{M}(m)");
2695   TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
2696
2697   TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
2698   // norm: 1/(sqrt(2pi)sigma)
2699   gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
2700   //gaus->Print();
2701
2702   for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
2703   {
2704     Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
2705     inputDist->SetBinContent(x, value);
2706     inputDistVector(x-1) = value;
2707   }
2708
2709   TVectorD measuredDistIdealVector = matrix * inputDistVector;
2710   
2711   for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
2712     measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
2713
2714   measuredDist->FillRandom(measuredIdealDist, 10000);
2715
2716   // fill error matrix before scaling
2717   TMatrixD covarianceMatrixMeasured(kSize, kSize);
2718   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
2719     covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
2720
2721   TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
2722   //covarianceMatrix.Print();
2723
2724   TVectorD measuredDistVector(kSize);
2725   for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
2726     measuredDistVector(x-1) = measuredDist->GetBinContent(x);
2727
2728   TVectorD unfoldedVector = inverted * measuredDistVector;
2729   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
2730     unfolded->SetBinContent(x, unfoldedVector(x-1));
2731
2732   TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1000, 500);
2733   canvas->SetTopMargin(0.05);
2734   canvas->Divide(2, 1);
2735
2736   canvas->cd(1);
2737   canvas->cd(1)->SetLeftMargin(0.15);
2738   canvas->cd(1)->SetRightMargin(0.05);
2739   measuredDist->GetYaxis()->SetTitleOffset(1.7);
2740   measuredDist->SetStats(0);
2741   measuredDist->DrawCopy();
2742   gaus->Draw("SAME");
2743
2744   canvas->cd(2);
2745   canvas->cd(2)->SetLeftMargin(0.15);
2746   canvas->cd(2)->SetRightMargin(0.05);
2747   unfolded->GetYaxis()->SetTitleOffset(1.7);
2748   unfolded->SetStats(0);
2749   unfolded->DrawCopy();
2750   gaus->Draw("SAME");
2751
2752   canvas->SaveAs("BlobelUnfoldingExample.eps");
2753 }
2754
2755 void E735Fit()
2756 {
2757   TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
2758   fCurrentESD->Sumw2();
2759
2760   // Open the input stream
2761   ifstream in;
2762   in.open("e735data.txt");
2763
2764   while(in.good())
2765   {
2766     Float_t x, y, ye;
2767     in >> x >> y >> ye;
2768
2769     //Printf("%f %f %f", x, y, ye);
2770     fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
2771     fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
2772   }
2773
2774   in.close();
2775
2776   //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
2777
2778   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2779
2780   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])");
2781   func->SetParNames("scaling", "averagen", "k");
2782   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
2783   func->SetParLimits(1, 0.001, 1000);
2784   func->SetParLimits(2, 0.001, 1000);
2785   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
2786
2787   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
2788   lognormal->SetParNames("scaling", "mean", "sigma");
2789   lognormal->SetParameters(1, 1, 1);
2790   lognormal->SetParLimits(0, 0, 10);
2791   lognormal->SetParLimits(1, 0, 100);
2792   lognormal->SetParLimits(2, 1e-3, 10);
2793
2794   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
2795   fCurrentESD->SetStats(kFALSE);
2796   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
2797   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
2798   fCurrentESD->Draw("");
2799   fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
2800   fCurrentESD->Fit(func, "0", "", 0, 150);
2801   func->SetRange(0, 250);
2802   func->Draw("SAME");
2803   printf("chi2 = %f\n", func->GetChisquare());
2804
2805   fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
2806   lognormal->SetLineColor(2);
2807   lognormal->SetLineStyle(2);
2808   lognormal->SetRange(0, 250);
2809   lognormal->Draw("SAME");
2810
2811   gPad->SetLogy();
2812
2813   canvas->SaveAs("E735Fit.eps");
2814 }