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