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