coding conventions
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / plotsMultiplicity.C
1 //
2 // plots for the note about multplicity measurements
3 //
4
5 #if !defined(__CINT__) || defined(__MAKECINT__)
6
7 #include <TCanvas.h>
8 #include <TPad.h>
9 #include <TH1F.h>
10 #include <TH2F.h>
11 #include <TH3F.h>
12 #include <TLine.h>
13 #include <TF1.h>
14 #include <TSystem.h>
15 #include <TFile.h>
16 #include <TLegend.h>
17 #include <TStopwatch.h>
18 #include <TROOT.h>
19 #include <TGraph.h>
20 #include <TMath.h>
21
22 #include "AliMultiplicityCorrection.h"
23 #include "AliCorrection.h"
24 #include "AliCorrectionMatrix3D.h"
25
26 #endif
27
28 const char* correctionFile = "multiplicityMC_2M.root";
29 const char* measuredFile   = "multiplicityMC_1M_3.root";
30 Int_t etaRange = 3;
31 Int_t drawRatioRange = 150;
32
33 const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
34 const char* measuredFileTPC   = "multiplicityMC_TPC_0.6M.root";
35 Int_t etaRangeTPC = 1;
36
37 void SetTPC()
38 {
39   correctionFile = correctionFileTPC;
40   measuredFile = measuredFileTPC;
41   etaRange = etaRangeTPC;
42   drawRatioRange = 100;
43 }
44
45 void Smooth(TH1* hist, Int_t windowWidth = 20)
46 {
47   TH1* clone = (TH1*) hist->Clone("clone");
48   for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
49   {
50     Int_t min = TMath::Max(2, bin-windowWidth);
51     Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
52     Float_t average = clone->Integral(min, max) / (max - min + 1);
53
54     hist->SetBinContent(bin, average);
55     hist->SetBinError(bin, 0);
56   }
57
58   delete clone;
59 }
60
61 void responseMatrixPlot()
62 {
63   gSystem->Load("libPWG0base");
64
65   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
66
67   TFile::Open(correctionFile);
68   mult->LoadHistograms("Multiplicity");
69
70   TH1* hist = mult->GetCorrelation(etaRange)->Project3D("zy");
71   hist->SetStats(kFALSE);
72
73   hist->SetTitle(";true multiplicity;measured multiplicity;Entries");
74   hist->GetXaxis()->SetRangeUser(0, 200);
75   hist->GetYaxis()->SetRangeUser(0, 200);
76
77   TCanvas* canvas = new TCanvas("c1", "c1", 800, 600);
78   canvas->SetRightMargin(0.15);
79   canvas->SetTopMargin(0.05);
80
81   gPad->SetLogz();
82   hist->Draw("COLZ");
83
84   canvas->SaveAs("responsematrix.eps");
85 }
86
87 TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
88 {
89   // normalize unfolded result to mc hist
90   result->Scale(1.0 / result->Integral(2, 200));
91   result->Scale(mcHist->Integral(2, 200));
92
93   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
94   canvas->Range(0, 0, 1, 1);
95
96   TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
97   pad1->Draw();
98
99   TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
100   pad2->Draw();
101
102   pad1->SetRightMargin(0.05);
103   pad2->SetRightMargin(0.05);
104
105   // no border between them
106   pad1->SetBottomMargin(0);
107   pad2->SetTopMargin(0);
108
109   pad1->cd();
110
111   mcHist->GetXaxis()->SetLabelSize(0.06);
112   mcHist->GetYaxis()->SetLabelSize(0.06);
113   mcHist->GetXaxis()->SetTitleSize(0.06);
114   mcHist->GetYaxis()->SetTitleSize(0.06);
115   mcHist->GetYaxis()->SetTitleOffset(0.6);
116
117   mcHist->GetXaxis()->SetRangeUser(0, 200);
118
119   mcHist->SetTitle(";true multiplicity;Entries");
120   mcHist->SetStats(kFALSE);
121
122   mcHist->DrawCopy("HIST E");
123   gPad->SetLogy();
124
125   result->SetLineColor(2);
126   result->DrawCopy("SAME HISTE");
127
128   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
129   legend->AddEntry(mcHist, "true distribution");
130   legend->AddEntry(result, "unfolded distribution");
131   legend->SetFillColor(0);
132   legend->Draw();
133
134   pad2->cd();
135   pad2->SetBottomMargin(0.15);
136
137   // calculate ratio
138   mcHist->Sumw2();
139   TH1* ratio = (TH1*) mcHist->Clone("ratio");
140   result->Sumw2();
141   ratio->Divide(ratio, result, 1, 1, "");
142   ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
143   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
144
145   ratio->DrawCopy();
146
147   // get average of ratio
148   Float_t sum = 0;
149   for (Int_t i=2; i<=150; ++i)
150   {
151     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
152   }
153   sum /= 149;
154
155   printf("Average (2..150) of |ratio - 1| is %f\n", sum);
156
157   TLine* line = new TLine(0, 1, 200, 1);
158   line->SetLineWidth(2);
159   line->Draw();
160
161   line = new TLine(0, 1.1, 200, 1.1);
162   line->SetLineWidth(2);
163   line->SetLineStyle(2);
164   line->Draw();
165   line = new TLine(0, 0.9, 200, 0.9);
166   line->SetLineWidth(2);
167   line->SetLineStyle(2);
168   line->Draw();
169
170   canvas->Modified();
171
172   canvas->SaveAs(epsName);
173
174   return canvas;
175 }
176
177 TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
178 {
179   // draws the 3 plots in the upper plot
180   // draws the ratio between result1 and result2 in the lower plot
181
182   // normalize unfolded result to mc hist
183   result1->Scale(1.0 / result1->Integral(2, 200));
184   result1->Scale(mcHist->Integral(2, 200));
185   result2->Scale(1.0 / result2->Integral(2, 200));
186   result2->Scale(mcHist->Integral(2, 200));
187
188   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
189   canvas->Range(0, 0, 1, 1);
190
191   TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
192   pad1->Draw();
193
194   TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
195   pad2->Draw();
196
197   pad1->SetRightMargin(0.05);
198   pad2->SetRightMargin(0.05);
199
200   // no border between them
201   pad1->SetBottomMargin(0);
202   pad2->SetTopMargin(0);
203
204   pad1->cd();
205
206   mcHist->GetXaxis()->SetLabelSize(0.06);
207   mcHist->GetYaxis()->SetLabelSize(0.06);
208   mcHist->GetXaxis()->SetTitleSize(0.06);
209   mcHist->GetYaxis()->SetTitleSize(0.06);
210   mcHist->GetYaxis()->SetTitleOffset(0.6);
211
212   mcHist->GetXaxis()->SetRangeUser(0, drawRatioRange);
213
214   mcHist->SetTitle(";true multiplicity;Entries");
215   mcHist->SetStats(kFALSE);
216
217   mcHist->DrawCopy("HIST E");
218   gPad->SetLogy();
219
220   result1->SetLineColor(2);
221   result1->DrawCopy("SAME HISTE");
222
223   result2->SetLineColor(4);
224   result2->DrawCopy("SAME HISTE");
225
226   TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
227   legend->AddEntry(mcHist, "true distribution");
228   legend->AddEntry(result1, "unfolded distribution (syst)");
229   legend->AddEntry(result2, "unfolded distribution (normal)");
230   legend->SetFillColor(0);
231   legend->Draw();
232
233   pad2->cd();
234   pad2->SetBottomMargin(0.15);
235
236   result1->GetXaxis()->SetLabelSize(0.06);
237   result1->GetYaxis()->SetLabelSize(0.06);
238   result1->GetXaxis()->SetTitleSize(0.06);
239   result1->GetYaxis()->SetTitleSize(0.06);
240   result1->GetYaxis()->SetTitleOffset(0.6);
241
242   result1->GetXaxis()->SetRangeUser(0, drawRatioRange);
243
244   result1->SetTitle(";true multiplicity;Entries");
245   result1->SetStats(kFALSE);
246
247   // calculate ratio
248   result1->Sumw2();
249   TH1* ratio = (TH1*) result1->Clone("ratio");
250   result2->Sumw2();
251   ratio->Divide(ratio, result2, 1, 1, "");
252   ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
253   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
254
255   ratio->DrawCopy();
256
257   // get average of ratio
258   Float_t sum = 0;
259   for (Int_t i=2; i<=drawRatioRange; ++i)
260   {
261     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
262   }
263   sum /= drawRatioRange-1;
264
265   printf("Average (2..%d) of |ratio - 1| is %f\n", drawRatioRange, sum);
266
267   TLine* line = new TLine(0, 1, drawRatioRange, 1);
268   line->SetLineWidth(2);
269   line->Draw();
270
271   line = new TLine(0, 1.1, drawRatioRange, 1.1);
272   line->SetLineWidth(2);
273   line->SetLineStyle(2);
274   line->Draw();
275   line = new TLine(0, 0.9, drawRatioRange, 0.9);
276   line->SetLineWidth(2);
277   line->SetLineStyle(2);
278   line->Draw();
279
280   canvas->Modified();
281
282   canvas->SaveAs(epsName);
283
284   return canvas;
285 }
286
287 TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0)
288 {
289   // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
290   // a systematic effect
291
292   // normalize results
293   result->Scale(1.0 / result->Integral(2, 200));
294
295   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
296
297   result->GetXaxis()->SetRangeUser(1, drawRatioRange);
298   result->SetStats(kFALSE);
299
300   Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
301
302   TLegend* legend = new TLegend(0.2, 0.75, 0.35, 0.95);
303   legend->SetFillColor(0);
304
305   for (Int_t n=0; n<nResultSyst; ++n)
306   {
307     resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
308
309     // calculate ratio
310     TH1* ratio = (TH1*) result->Clone("ratio");
311     ratio->Divide(ratio, resultSyst[n], 1, 1, "B");
312     ratio->SetTitle(";true multiplicity;Ratio");
313     ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
314
315     if (firstMarker)
316       ratio->SetMarkerStyle(5);
317
318     ratio->SetLineColor(colors[n / 2]);
319     if ((n % 2))
320       ratio->SetLineStyle(2);
321
322     ratio->DrawCopy((n == 0) ? ((firstMarker) ? "P" : "HIST") : "SAME HIST");
323
324     if (legendStrings && legendStrings[n])
325       legend->AddEntry(ratio, legendStrings[n]);
326
327     // get average of ratio
328     Float_t sum = 0;
329     for (Int_t i=2; i<=drawRatioRange; ++i)
330       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
331     sum /= drawRatioRange-1;
332
333     printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, drawRatioRange, sum);
334   }
335
336   if (legendStrings)
337     legend->Draw();
338
339   TLine* line = new TLine(1, 1, drawRatioRange, 1);
340   line->SetLineWidth(2);
341   line->Draw();
342
343   line = new TLine(1, 1.1, drawRatioRange, 1.1);
344   line->SetLineWidth(2);
345   line->SetLineStyle(2);
346   line->Draw();
347   line = new TLine(1, 0.9, drawRatioRange, 0.9);
348   line->SetLineWidth(2);
349   line->SetLineStyle(2);
350   line->Draw();
351
352   canvas->Modified();
353
354   canvas->SaveAs(epsName);
355   canvas->SaveAs(Form("%s.gif", epsName.Data()));
356
357   return canvas;
358 }
359
360 TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
361 {
362   // draws the ratios of each mc to the corresponding result
363
364   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
365   canvas->SetRightMargin(0.05);
366   canvas->SetTopMargin(0.05);
367
368   for (Int_t n=0; n<nResultSyst; ++n)
369   {
370     // normalize
371     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
372     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
373
374     result[n]->GetXaxis()->SetRangeUser(0, drawRatioRange);
375     result[n]->SetStats(kFALSE);
376
377     // calculate ratio
378     TH1* ratio = (TH1*) result[n]->Clone("ratio");
379     ratio->Divide(mc[n], ratio, 1, 1, "B");
380
381     // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
382     ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
383
384     if (smooth)
385       Smooth(ratio);
386
387     ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
388     ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
389
390     if (dashed)
391     {
392       ratio->SetLineColor((n/2)+1);
393       ratio->SetLineStyle((n%2)+1);
394     }
395     else
396       ratio->SetLineColor(n+1);
397
398     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
399
400     // get average of ratio
401     Float_t sum = 0;
402     for (Int_t i=2; i<=drawRatioRange; ++i)
403       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
404     sum /= drawRatioRange-1;
405
406     printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, drawRatioRange, sum);
407   }
408
409   TLine* line = new TLine(0, 1, drawRatioRange, 1);
410   line->SetLineWidth(2);
411   line->Draw();
412
413   line = new TLine(0, 1.1, drawRatioRange, 1.1);
414   line->SetLineWidth(2);
415   line->SetLineStyle(2);
416   line->Draw();
417   line = new TLine(0, 0.9, drawRatioRange, 0.9);
418   line->SetLineWidth(2);
419   line->SetLineStyle(2);
420   line->Draw();
421
422   canvas->Modified();
423
424   canvas->SaveAs(epsName);
425   canvas->SaveAs(Form("%s.gif", epsName.Data()));
426
427   return canvas;
428 }
429
430 TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
431 {
432   // draws the ratios of each mc to the corresponding result
433   // deducts from each ratio the ratio of mcBase / resultBase
434
435   // normalize
436   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
437   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
438
439   // calculate ratio
440   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
441   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
442
443   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
444   canvas->SetRightMargin(0.05);
445   canvas->SetTopMargin(0.05);
446
447   for (Int_t n=0; n<nResultSyst; ++n)
448   {
449     // normalize
450     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
451     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
452
453     result[n]->GetXaxis()->SetRangeUser(0, drawRatioRange);
454     result[n]->SetStats(kFALSE);
455
456     // calculate ratio
457     TH1* ratio = (TH1*) result[n]->Clone("ratio");
458     ratio->Divide(mc[n], ratio, 1, 1, "B");
459     ratio->Add(ratioBase, -1);
460
461     ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
462     ratio->GetYaxis()->SetRangeUser(-1, 1);
463     ratio->SetLineColor(n+1);
464     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
465
466     // get average of ratio
467     Float_t sum = 0;
468     for (Int_t i=2; i<=drawRatioRange; ++i)
469       sum += TMath::Abs(ratio->GetBinContent(i));
470     sum /= drawRatioRange-1;
471
472     printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, drawRatioRange, sum);
473   }
474
475   TLine* line = new TLine(0, 0, drawRatioRange, 0);
476   line->SetLineWidth(2);
477   line->Draw();
478
479   line = new TLine(0, 0.1, drawRatioRange, 0.1);
480   line->SetLineWidth(2);
481   line->SetLineStyle(2);
482   line->Draw();
483   line = new TLine(0, -0.1, drawRatioRange, -0.1);
484   line->SetLineWidth(2);
485   line->SetLineStyle(2);
486   line->Draw();
487
488   canvas->Modified();
489
490   canvas->SaveAs(epsName);
491   canvas->SaveAs(Form("%s.gif", epsName.Data()));
492
493   return canvas;
494 }
495
496 TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
497 {
498   // draws the ratios of each mc to the corresponding result
499   // deducts from each ratio the ratio of mcBase / resultBase
500   // smoothens the ratios by a sliding window
501
502   // normalize
503   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
504   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
505
506   // calculate ratio
507   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
508   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
509
510   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
511   canvas->SetRightMargin(0.05);
512   canvas->SetTopMargin(0.05);
513
514   for (Int_t n=0; n<nResultSyst; ++n)
515   {
516     // normalize
517     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
518     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
519
520     result[n]->GetXaxis()->SetRangeUser(0, drawRatioRange);
521     result[n]->SetStats(kFALSE);
522
523     // calculate ratio
524     TH1* ratio = (TH1*) result[n]->Clone("ratio");
525     ratio->Divide(mc[n], ratio, 1, 1, "B");
526     ratio->Add(ratioBase, -1);
527
528     //new TCanvas; ratio->DrawCopy();
529     // clear 0 bin
530     ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
531
532     Smooth(ratio);
533
534     //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
535
536     canvas->cd();
537     ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
538     ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
539     ratio->SetLineColor((n / 2)+1);
540     ratio->SetLineStyle((n % 2)+1);
541     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
542
543     // get average of ratio
544     Float_t sum = 0;
545     for (Int_t i=2; i<=150; ++i)
546       sum += TMath::Abs(ratio->GetBinContent(i));
547     sum /= 149;
548
549     printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
550   }
551
552   TLine* line = new TLine(0, 0, drawRatioRange, 0);
553   line->SetLineWidth(2);
554   line->Draw();
555
556   line = new TLine(0, 0.1, drawRatioRange, 0.1);
557   line->SetLineWidth(2);
558   line->SetLineStyle(2);
559   line->Draw();
560   line = new TLine(0, -0.1, drawRatioRange, -0.1);
561   line->SetLineWidth(2);
562   line->SetLineStyle(2);
563   line->Draw();
564
565   canvas->Modified();
566
567   canvas->SaveAs(epsName);
568   canvas->SaveAs(Form("%s.gif", epsName.Data()));
569
570   return canvas;
571 }
572
573 void DrawResiduals(TH1* measured, TH1* unfoldedFolded, const char* epsName)
574 {
575   // normalize
576   unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, 200));
577   unfoldedFolded->Scale(measured->Integral(2, 200));
578
579   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
580   canvas->Range(0, 0, 1, 1);
581
582   TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
583   pad1->Draw();
584   pad1->SetGridx();
585   pad1->SetGridy();
586
587   TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
588   pad2->Draw();
589   pad2->SetGridx();
590   pad2->SetGridy();
591
592   TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
593   pad3->SetGridx();
594   pad3->SetGridy();
595   pad3->SetRightMargin(0.05);
596   pad3->SetTopMargin(0.05);
597   pad3->Draw();
598
599   pad1->SetRightMargin(0.05);
600   pad2->SetRightMargin(0.05);
601
602   // no border between them
603   pad1->SetBottomMargin(0);
604   pad2->SetTopMargin(0);
605
606   pad1->cd();
607
608   measured->GetXaxis()->SetLabelSize(0.06);
609   measured->GetYaxis()->SetLabelSize(0.06);
610   measured->GetXaxis()->SetTitleSize(0.06);
611   measured->GetYaxis()->SetTitleSize(0.06);
612   measured->GetYaxis()->SetTitleOffset(0.6);
613
614   measured->GetXaxis()->SetRangeUser(0, 150);
615
616   measured->SetTitle(";measured multiplicity;Entries");
617   measured->SetStats(kFALSE);
618
619   measured->DrawCopy("HIST");
620   gPad->SetLogy();
621
622   unfoldedFolded->SetMarkerStyle(5);
623   unfoldedFolded->SetMarkerColor(2);
624   unfoldedFolded->SetLineColor(0);
625   unfoldedFolded->DrawCopy("SAME P");
626
627   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
628   legend->AddEntry(measured, "measured distribution");
629   legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution");
630   legend->SetFillColor(0);
631   legend->Draw();
632
633   pad2->cd();
634   pad2->SetBottomMargin(0.15);
635
636   // calculate ratio
637   measured->Sumw2();
638   TH1* residual = (TH1*) measured->Clone("residual");
639   unfoldedFolded->Sumw2();
640
641   residual->Add(unfoldedFolded, -1);
642
643   // projection
644   TH1* residualHist = new TH1F("residualHist", ";", 15, -3, 3);
645
646   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
647   {
648     if (measured->GetBinError(i) > 0)
649     {
650       residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
651       residual->SetBinError(i, 1);
652
653       residualHist->Fill(residual->GetBinContent(i));
654     }
655     else
656     {
657       residual->SetBinContent(i, 0);
658       residual->SetBinError(i, 0);
659     }
660   }
661
662   residual->GetYaxis()->SetTitle("Residuals   1/e (M - R #otimes U)");
663   residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
664   residual->DrawCopy();
665
666   TLine* line = new TLine(-0.5, 0, 150.5, 0);
667   line->SetLineWidth(2);
668   line->Draw();
669
670   pad3->cd();
671   residualHist->SetStats(kFALSE);
672   residualHist->GetXaxis()->SetLabelSize(0.08);
673   residualHist->Fit("gaus");
674   residualHist->Draw();
675
676   canvas->Modified();
677   canvas->SaveAs(canvas->GetName());
678
679   //const char* epsName2 = "proj.eps";
680   //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
681   //canvas->SetGridx();
682   //canvas->SetGridy();
683
684   //canvas->SaveAs(canvas->GetName());
685 }
686
687 void bayesianExample()
688 {
689   TStopwatch watch;
690   watch.Start();
691
692   gSystem->Load("libPWG0base");
693
694   TFile::Open(correctionFile);
695   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
696   mult->LoadHistograms("Multiplicity");
697
698   TFile::Open(measuredFile);
699   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
700   mult2->LoadHistograms("Multiplicity");
701
702   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
703
704   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
705
706   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
707   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
708
709   //mult->DrawComparison("bayesianExample", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
710   DrawResultRatio(mcHist, result, "bayesianExample.eps");
711
712   // draw residual plot
713
714   // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
715   TH2* convoluted = mult->CalculateMultiplicityESD(result, etaRange);
716   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
717
718   TH1* measured = mult2->GetMultiplicityESD(etaRange)->ProjectionY("measured");
719
720   DrawResiduals(measured, convolutedProj, "bayesianResiduals.eps");
721
722   watch.Stop();
723   watch.Print();
724 }
725
726 void chi2FluctuationResult()
727 {
728   gSystem->Load("libPWG0base");
729
730   TFile::Open(correctionFile);
731   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
732   mult->LoadHistograms("Multiplicity");
733
734   TFile::Open(measuredFile);
735   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
736   mult2->LoadHistograms("Multiplicity");
737
738   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
739   mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
740   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
741
742   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
743   //TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
744
745   mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
746
747   TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_3");
748   canvas->SaveAs("chi2FluctuationResult.eps");
749 }
750
751 void chi2Example()
752 {
753   TStopwatch watch;
754   watch.Start();
755
756   gSystem->Load("libPWG0base");
757
758   TFile::Open(correctionFile);
759   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
760   mult->LoadHistograms("Multiplicity");
761
762   TFile::Open(measuredFile);
763   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
764   mult2->LoadHistograms("Multiplicity");
765
766   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
767   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
768   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
769
770   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
771   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
772
773   DrawResultRatio(mcHist, result, "chi2Example.eps");
774
775   watch.Stop();
776   watch.Print();
777 }
778
779 void chi2ExampleTPC()
780 {
781   TStopwatch watch;
782   watch.Start();
783
784   gSystem->Load("libPWG0base");
785
786   TFile::Open(correctionFileTPC);
787   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
788   mult->LoadHistograms("Multiplicity");
789
790   TFile::Open(measuredFileTPC);
791   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
792   mult2->LoadHistograms("Multiplicity");
793
794   mult->SetMultiplicityESD(etaRangeTPC, mult2->GetMultiplicityESD(etaRangeTPC));
795   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
796   mult->ApplyMinuitFit(etaRangeTPC, kFALSE, AliMultiplicityCorrection::kTrVtx);
797
798   TH1* mcHist = mult2->GetMultiplicityVtx(etaRangeTPC)->ProjectionY("mymc");
799   TH1* result = mult->GetMultiplicityESDCorrected(etaRangeTPC);
800
801   DrawResultRatio(mcHist, result, "chi2ExampleTPC.eps");
802
803   watch.Stop();
804   watch.Print();
805 }
806
807 void bayesianNBD()
808 {
809   gSystem->Load("libPWG0base");
810   TFile::Open("multiplicityMC_3M.root");
811
812   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
813   mult->LoadHistograms("Multiplicity");
814
815   TFile::Open("multiplicityMC_3M_NBD.root");
816   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
817   mult2->LoadHistograms("Multiplicity");
818
819   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
820   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1.0, 100);
821
822   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
823
824   mcHist->Sumw2();
825   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
826
827   //mult->DrawComparison("bayesianNBD", etaRange, kFALSE, kTRUE, mcHist);
828   DrawResultRatio(mcHist, result, "bayesianNBD.eps");
829 }
830
831 void minimizationNBD()
832 {
833   gSystem->Load("libPWG0base");
834   TFile::Open("multiplicityMC_3M.root");
835
836   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
837   mult->LoadHistograms("Multiplicity");
838
839   TFile::Open("multiplicityMC_3M_NBD.root");
840   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
841   mult2->LoadHistograms("Multiplicity");
842
843   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
844   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
845   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
846
847   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
848
849   mcHist->Sumw2();
850   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
851
852   //mult->DrawComparison("minimizationNBD", etaRange, kFALSE, kTRUE, mcHist);
853   DrawResultRatio(mcHist, result, "minimizationNBD.eps");
854 }
855
856 void minimizationInfluenceAlpha()
857 {
858   gSystem->Load("libPWG0base");
859
860   TFile::Open(measuredFile);
861   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
862   mult2->LoadHistograms("Multiplicity");
863
864   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
865   mcHist->Scale(1.0 / mcHist->Integral());
866   mcHist->GetXaxis()->SetRangeUser(0, 200);
867   mcHist->SetStats(kFALSE);
868   mcHist->SetTitle(";true multiplicity;P_{N}");
869
870   TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 300);
871   canvas->Divide(3, 1);
872
873   TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
874
875   TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_100.000000");
876   TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_03_2_100000.000000");
877   TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_06_2_100000000.000000");
878
879   mcHist->Rebin(2);  mcHist->Scale(0.5);
880   hist1->Rebin(2);   hist1->Scale(0.5);
881   hist2->Rebin(2);   hist2->Scale(0.5);
882   hist3->Rebin(2);   hist3->Scale(0.5);
883
884   mcHist->GetXaxis()->SetRangeUser(0, 200);
885
886   canvas->cd(1);
887   gPad->SetLogy();
888   mcHist->Draw();
889   hist1->SetMarkerStyle(5);
890   hist1->SetMarkerColor(2);
891   hist1->Draw("SAME PE");
892
893   canvas->cd(2);
894   gPad->SetLogy();
895   mcHist->Draw();
896   hist2->SetMarkerStyle(5);
897   hist2->SetMarkerColor(2);
898   hist2->Draw("SAME PE");
899
900   canvas->cd(3);
901   gPad->SetLogy();
902   mcHist->Draw();
903   hist3->SetMarkerStyle(5);
904   hist3->SetMarkerColor(2);
905   hist3->Draw("SAME PE");
906
907   canvas->SaveAs("minimizationInfluenceAlpha.eps");
908 }
909
910 void NBDFit()
911 {
912   gSystem->Load("libPWG0base");
913
914   TFile::Open(correctionFile);
915   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
916   mult->LoadHistograms("Multiplicity");
917
918   TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
919   fCurrentESD->Sumw2();
920   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
921
922   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])");
923   func->SetParNames("scaling", "averagen", "k");
924   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
925   func->SetParLimits(1, 0.001, 1000);
926   func->SetParLimits(2, 0.001, 1000);
927   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
928
929   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
930   lognormal->SetParNames("scaling", "mean", "sigma");
931   lognormal->SetParameters(1, 1, 1);
932   lognormal->SetParLimits(0, 0, 10);
933   lognormal->SetParLimits(1, 0, 100);
934   lognormal->SetParLimits(2, 1e-3, 10);
935
936   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
937   fCurrentESD->SetStats(kFALSE);
938   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
939   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
940   fCurrentESD->Draw("HIST");
941   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
942   fCurrentESD->Fit(func, "W0", "", 0, 50);
943   func->SetRange(0, 100);
944   func->Draw("SAME");
945   printf("chi2 = %f\n", func->GetChisquare());
946
947   fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
948   lognormal->SetLineColor(2);
949   lognormal->SetLineStyle(2);
950   lognormal->SetRange(0, 100);
951   lognormal->Draw("SAME");
952
953   canvas->SaveAs("NBDFit.eps");
954 }
955
956 void DifferentSamples()
957 {
958   // data generated by runMultiplicitySelector.C DifferentSamples
959
960   const char* name = "DifferentSamples";
961
962   TFile* file = TFile::Open(Form("%s.root", name));
963
964   TCanvas* canvas = new TCanvas(name, name, 800, 600);
965   canvas->Divide(2, 2);
966
967   for (Int_t i=0; i<4; ++i)
968   {
969     canvas->cd(i+1);
970     gPad->SetTopMargin(0.05);
971     gPad->SetRightMargin(0.05);
972     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", i));
973     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", i));
974     TH1* mc = (TH1*) file->Get(Form("mc_%d", i));
975     mc->Sumw2();
976
977     chi2Result->Divide(chi2Result, mc, 1, 1, "");
978     bayesResult->Divide(bayesResult, mc, 1, 1, "");
979
980     chi2Result->SetTitle(";true multiplicity;unfolded measured/MC");
981     chi2Result->GetXaxis()->SetRangeUser(0, 150);
982     chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
983     chi2Result->GetYaxis()->SetTitleOffset(1.2);
984     chi2Result->SetLineColor(1);
985     chi2Result->SetStats(kFALSE);
986
987     bayesResult->SetStats(kFALSE);
988     bayesResult->SetLineColor(2);
989
990     chi2Result->DrawCopy("HIST");
991     bayesResult->DrawCopy("SAME HIST");
992
993     TLine* line = new TLine(0, 1, 150, 1);
994     line->Draw();
995   }
996
997   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
998 }
999
1000 void StartingConditions()
1001 {
1002   // data generated by runMultiplicitySelector.C StartingConditions
1003
1004   const char* name = "StartingConditions";
1005
1006   TFile* file = TFile::Open(Form("%s.root", name));
1007
1008   TCanvas* canvas = new TCanvas(name, name, 800, 400);
1009   canvas->Divide(2, 1);
1010
1011   TH1* mc = (TH1*) file->Get("mc");
1012   mc->Sumw2();
1013   mc->Scale(1.0 / mc->Integral());
1014
1015   //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
1016
1017   for (Int_t i=0; i<6; ++i)
1018   {
1019     Int_t id = i;
1020     if (id > 2)
1021       id += 2;
1022
1023     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1024     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
1025
1026     chi2Result->Divide(chi2Result, mc, 1, 1, "");
1027     bayesResult->Divide(bayesResult, mc, 1, 1, "");
1028
1029     chi2Result->SetTitle("a) #chi^{2} minimization;true multiplicity;unfolded / MC");
1030     chi2Result->GetXaxis()->SetRangeUser(0, 150);
1031     chi2Result->GetYaxis()->SetRangeUser(0.8, 1.2);
1032     chi2Result->GetYaxis()->SetTitleOffset(1.5);
1033     //chi2Result->SetMarkerStyle(marker[i]);
1034     chi2Result->SetLineColor(i+1);
1035     chi2Result->SetStats(kFALSE);
1036
1037     bayesResult->SetTitle("b) Bayesian method;true multiplicity;unfolded / MC");
1038     bayesResult->GetXaxis()->SetRangeUser(0, 150);
1039     bayesResult->GetYaxis()->SetRangeUser(0.8, 1.2);
1040     bayesResult->GetYaxis()->SetTitleOffset(1.5);
1041     bayesResult->SetStats(kFALSE);
1042     bayesResult->SetLineColor(2);
1043     bayesResult->SetLineColor(i+1);
1044
1045     canvas->cd(1);
1046     gPad->SetLeftMargin(0.12);
1047     chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1048
1049     canvas->cd(2);
1050     gPad->SetLeftMargin(0.12);
1051     bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1052
1053     //TLine* line = new TLine(0, 1, 150, 1);
1054     //line->Draw();
1055   }
1056
1057   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1058 }
1059
1060 void StatisticsPlot()
1061 {
1062   const char* name = "StatisticsPlot";
1063
1064   TFile* file = TFile::Open(Form("%s.root", name));
1065
1066   TCanvas* canvas = new TCanvas(name, name, 600, 400);
1067
1068   TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
1069   fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1070   fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1071   fitResultsChi2->Draw("AP");
1072
1073   TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1074   fitResultsChi2->Fit(f);
1075
1076   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1077
1078   TH1* mc[5];
1079   TH1* result[5];
1080
1081   const char* plotname = "chi2Result";
1082
1083   name = "StatisticsPlotRatios";
1084   canvas = new TCanvas(name, name, 600, 400);
1085
1086   for (Int_t i=0; i<5; ++i)
1087   {
1088     mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1089     result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1090
1091     result[i]->SetLineColor(i+1);
1092     result[i]->Draw(((i == 0) ? "" : "SAME"));
1093   }
1094 }
1095
1096 void SystematicLowEfficiency()
1097 {
1098   gSystem->Load("libPWG0base");
1099
1100   TFile::Open(correctionFile);
1101   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1102   mult->LoadHistograms("Multiplicity");
1103
1104   // calculate result with systematic effect
1105   TFile::Open("multiplicityMC_100k_1_loweff.root");
1106   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1107   mult2->LoadHistograms("Multiplicity");
1108
1109   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1110   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1111   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1112
1113   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1114   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1115
1116   DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
1117
1118   // calculate normal result
1119   TFile::Open("multiplicityMC_100k_1.root");
1120   mult2->LoadHistograms("Multiplicity");
1121
1122   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1123   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1124
1125   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1126
1127   DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
1128
1129   Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
1130 }
1131
1132 void SystematicMisalignment()
1133 {
1134   gSystem->Load("libPWG0base");
1135
1136   TFile::Open(correctionFile);
1137   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1138   mult->LoadHistograms("Multiplicity");
1139
1140   TFile::Open("multiplicityMC_100k_fullmis.root");
1141   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1142   mult2->LoadHistograms("Multiplicity");
1143
1144   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1145   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1146   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1147
1148   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1149   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1150
1151   DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1152 }
1153
1154 void SystematicMisalignmentTPC()
1155 {
1156   gSystem->Load("libPWG0base");
1157
1158   SetTPC();
1159
1160   TFile::Open(correctionFile);
1161   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1162   mult->LoadHistograms("Multiplicity");
1163
1164   TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1165   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1166   mult2->LoadHistograms("Multiplicity");
1167
1168   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1169   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1170   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1171
1172   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1173   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1174
1175   DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1176 }
1177
1178 void EfficiencySpecies()
1179 {
1180   gSystem->Load("libPWG0base");
1181
1182   Int_t marker[] = {24, 25, 26};
1183   Int_t color[] = {1, 2, 4};
1184
1185   // SPD TPC
1186   const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
1187   Float_t etaRange[] = {0.49, 0.9};
1188   const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1189
1190   TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1191   canvas->Divide(2, 1);
1192
1193   for (Int_t loop=0; loop<2; ++loop)
1194   {
1195     Printf("%s", fileName[loop]);
1196
1197     TFile::Open(fileName[loop]);
1198     AliCorrection* correction[4];
1199
1200     canvas->cd(loop+1);
1201
1202     gPad->SetGridx();
1203     gPad->SetGridy();
1204     gPad->SetRightMargin(0.05);
1205     //gPad->SetTopMargin(0.05);
1206
1207     TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
1208     legend->SetFillColor(0);
1209     legend->SetEntrySeparation(0.2);
1210
1211     Float_t below = 0;
1212     Float_t total = 0;
1213
1214     for (Int_t i=0; i<3; ++i)
1215     {
1216       Printf("correction %d", i);
1217
1218       TString name; name.Form("correction_%d", i);
1219       correction[i] = new AliCorrection(name, name);
1220       correction[i]->LoadHistograms();
1221
1222       TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1223       TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1224
1225       // limit vtx axis
1226       gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
1227       meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
1228
1229       // 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)
1230       /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1231         for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1232         {
1233           gene->SetBinContent(x, 0, z, 0);
1234           gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1235           meas->SetBinContent(x, 0, z, 0);
1236           meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1237         }*/
1238
1239       // limit eta axis
1240       gene->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1241       meas->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1242
1243       TH1* genePt = gene->Project3D(Form("z_%d", i));
1244       TH1* measPt = meas->Project3D(Form("z_%d", i));
1245
1246       genePt->Sumw2();
1247       measPt->Sumw2();
1248
1249       TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1250       effPt->Reset();
1251       effPt->Divide(measPt, genePt, 1, 1, "B");
1252
1253       Int_t bin = 0;
1254       for (bin=20; bin>=1; bin--)
1255       {
1256         if (effPt->GetBinContent(bin) < 0.5)
1257           break;
1258       }
1259
1260       Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1261
1262       Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1263       Printf("%.4f of the particles are below that momentum", fraction);
1264
1265       below += genePt->Integral(1, bin);
1266       total += genePt->Integral();
1267
1268       effPt->SetLineColor(color[i]);
1269       effPt->SetMarkerColor(color[i]);
1270       effPt->SetMarkerStyle(marker[i]);
1271
1272       effPt->GetXaxis()->SetRangeUser(0.06, 1);
1273       effPt->GetYaxis()->SetRangeUser(0, 1);
1274
1275       effPt->GetYaxis()->SetTitleOffset(1.2);
1276
1277       effPt->SetStats(kFALSE);
1278       effPt->SetTitle(titles[loop]);
1279       effPt->GetYaxis()->SetTitle("Efficiency");
1280
1281       effPt->DrawCopy((i == 0) ? "" : "SAME");
1282
1283       legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")));
1284     }
1285
1286     Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1287
1288     legend->Draw();
1289   }
1290
1291   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1292 }
1293
1294 void ParticleSpeciesComparison1(const char* fileNameMC = "multiplicityMC_400k_syst_species.root", const char* fileNameESD = "multiplicityMC_100k_syst.root")
1295 {
1296   gSystem->Load("libPWG0base");
1297
1298   Bool_t chi2 = 1;
1299
1300   TFile::Open(fileNameESD);
1301   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
1302   TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
1303
1304   TH1* results[10];
1305
1306   // loop over cases (normal, enhanced/reduced ratios)
1307   Int_t nMax = 7;
1308   for (Int_t i = 0; i<nMax; ++i)
1309   {
1310     TString folder;
1311     folder.Form("Multiplicity_%d", i);
1312
1313     AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
1314
1315     TFile::Open(fileNameMC);
1316     mult->LoadHistograms();
1317
1318     mult->SetMultiplicityESD(etaRange, hist);
1319
1320     if (chi2)
1321     {
1322       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1323       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1324       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1325     }
1326     else
1327     {
1328       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1329       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1330     }
1331
1332     //Float_t averageRatio = 0;
1333     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1334
1335     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1336
1337     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1338   }
1339
1340   DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
1341
1342   TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
1343
1344   for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1345   {
1346     results[0]->SetBinError(i, 0);
1347     mc->SetBinError(i, 0);
1348   }
1349
1350   const char* legendStrings[] = { "#pi^{#pm}", 0, "K^{#pm}", 0, "p,#bar{p}", 0 };
1351
1352   DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1353
1354   //not valid: draw chi2 uncertainty on top!
1355   /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1356   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1357   errorHist->SetLineColor(1);
1358   errorHist->SetLineWidth(2);
1359   TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1360   for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1361   {
1362     errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1363     errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1364   }
1365
1366   errorHist->DrawCopy("SAME");
1367   errorHist2->DrawCopy("SAME");*/
1368
1369   //canvas->SaveAs(canvas->GetName());
1370
1371   DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1372
1373   //errorHist->DrawCopy("SAME");
1374   //errorHist2->DrawCopy("SAME");
1375
1376   //canvas2->SaveAs(canvas2->GetName());
1377 }
1378
1379 void ParticleSpeciesComparison2()
1380 {
1381   gSystem->Load("libPWG0base");
1382
1383   const char* fileNameMC = "multiplicityMC_400k_syst.root";
1384   const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1385   Bool_t chi2 = 0;
1386
1387   TFile::Open(fileNameMC);
1388   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1389   mult->LoadHistograms();
1390
1391   TH1* mc[10];
1392   TH1* results[10];
1393
1394   // loop over cases (normal, enhanced/reduced ratios)
1395   Int_t nMax = 7;
1396   for (Int_t i = 0; i<nMax; ++i)
1397   {
1398     TString folder;
1399     folder.Form("Multiplicity_%d", i);
1400
1401     AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1402
1403     TFile::Open(fileNameESD);
1404     mult2->LoadHistograms();
1405
1406     mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1407
1408     if (chi2)
1409     {
1410       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1411       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1412       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1413     }
1414     else
1415     {
1416       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1417       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1418     }
1419
1420     //Float_t averageRatio = 0;
1421     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1422
1423     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1424
1425     TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1426     mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1427
1428     //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1429     //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1430
1431     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1432   }
1433
1434   DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1435 }
1436
1437 TH1* Invert(TH1* eff)
1438 {
1439   // calculate corr = 1 / eff
1440
1441   TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1442   corr->Reset();
1443
1444   for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1445   {
1446     if (eff->GetBinContent(i) > 0)
1447     {
1448       corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1449       corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1450     }
1451   }
1452
1453   return corr;
1454 }
1455
1456 void TriggerVertexCorrection()
1457 {
1458   //
1459   // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1460   //
1461
1462   gSystem->Load("libPWG0base");
1463
1464   TFile::Open(correctionFile);
1465   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1466   mult->LoadHistograms("Multiplicity");
1467
1468   TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1469   TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1470
1471   TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
1472
1473   corrINEL->SetStats(kFALSE);
1474   corrINEL->GetXaxis()->SetRangeUser(0, 30);
1475   corrINEL->GetYaxis()->SetRangeUser(0, 5);
1476   corrINEL->SetTitle(";true multiplicity;correction factor");
1477   corrINEL->SetMarkerStyle(22);
1478   corrINEL->Draw("PE");
1479
1480   corrMB->SetStats(kFALSE);
1481   corrMB->SetLineColor(2);
1482   corrMB->SetMarkerStyle(25);
1483   corrMB->SetMarkerColor(2);
1484   corrMB->Draw("SAME PE");
1485
1486   TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
1487   legend->SetFillColor(0);
1488   legend->AddEntry(corrINEL, "correction to inelastic sample");
1489   legend->AddEntry(corrMB, "correction to minimum bias sample");
1490
1491   legend->Draw();
1492
1493   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1494 }
1495
1496 void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1497 {
1498   gSystem->Load("libPWG0base");
1499
1500   TFile::Open(correctionFile);
1501   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1502   mult->LoadHistograms("Multiplicity");
1503
1504   TFile::Open(measuredFile);
1505   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1506   mult2->LoadHistograms("Multiplicity");
1507
1508   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1509
1510   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1511
1512   TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1513
1514   TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1515   TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1516
1517   if (!mc)
1518   {
1519     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1520     DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1521   }
1522
1523   TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1524   canvas->SetGridx();
1525   canvas->SetGridy();
1526   canvas->SetRightMargin(0.05);
1527   canvas->SetTopMargin(0.05);
1528
1529   errorResponse->SetLineColor(1);
1530   errorResponse->GetXaxis()->SetRangeUser(0, 200);
1531   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1532   errorResponse->SetStats(kFALSE);
1533   errorResponse->SetTitle(";true multiplicity;Uncertainty");
1534
1535   errorResponse->Draw();
1536
1537   errorMeasured->SetLineColor(2);
1538   errorMeasured->Draw("SAME");
1539
1540   errorBoth->SetLineColor(4);
1541   errorBoth->Draw("SAME");
1542
1543   Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
1544   Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
1545   Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
1546
1547   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1548
1549   TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
1550   errorResponse->Write();
1551   errorMeasured->Write();
1552   errorBoth->Write();
1553   file->Close();
1554 }
1555
1556 void StatisticalUncertaintyCompare(const char* det = "SPD")
1557 {
1558   TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1559   TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1560   TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1561   TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1562
1563   TString str;
1564   str.Form("StatisticalUncertaintyCompare%s", det);
1565
1566   TCanvas* canvas = new TCanvas(str, str, 600, 400);
1567   canvas->SetGridx();
1568   canvas->SetGridy();
1569   canvas->SetRightMargin(0.05);
1570   canvas->SetTopMargin(0.05);
1571
1572   errorResponse->SetLineColor(1);
1573   errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
1574   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1575   errorResponse->SetStats(kFALSE);
1576   errorResponse->SetTitle(";true multiplicity;Uncertainty");
1577
1578   errorResponse->Draw();
1579
1580   errorMeasured->SetLineColor(2);
1581   errorMeasured->Draw("SAME");
1582
1583   errorBoth->SetLineColor(4);
1584   errorBoth->Draw("SAME");
1585
1586   TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
1587   TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
1588
1589   errorBoth2->SetLineColor(4);
1590   errorBoth2->SetLineStyle(2);
1591   errorBoth2->Draw("SAME");
1592
1593   TLegend* legend = new TLegend(0.2, 0.6, 0.6, 0.9);
1594   legend->SetFillColor(0);
1595   legend->AddEntry(errorResponse, "response matrix (Bayesian)");
1596   legend->AddEntry(errorMeasured, "measured (Bayesian)");
1597   legend->AddEntry(errorBoth, "both (Bayesian)");
1598   legend->AddEntry(errorBoth2, "both (#chi^{2} minimization)");
1599   legend->Draw();
1600
1601   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1602 }
1603
1604 void EfficiencyComparison(Int_t eventType = 2)
1605 {
1606  const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
1607
1608   gSystem->Load("libPWG0base");
1609
1610   TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
1611   canvas->SetGridx();
1612   canvas->SetGridy();
1613   canvas->SetRightMargin(0.05);
1614   canvas->SetTopMargin(0.05);
1615
1616   AliMultiplicityCorrection* data[4];
1617   TH1* effArray[4];
1618
1619   Int_t markers[] = { 24, 25, 26, 5 };
1620   Int_t colors[] = { 1, 2, 3, 4 };
1621
1622   TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
1623   legend->SetFillColor(0);
1624
1625   TH1* effError = 0;
1626
1627   for (Int_t i=0; i<4; ++i)
1628   {
1629     TString name;
1630     name.Form("Multiplicity_%d", i);
1631
1632     TFile::Open(files[i]);
1633     data[i] = new AliMultiplicityCorrection(name, name);
1634
1635     if (i < 3)
1636     {
1637       data[i]->LoadHistograms("Multiplicity");
1638     }
1639     else
1640       data[i]->LoadHistograms("Multiplicity_0");
1641
1642     TH1* eff = (TH1*) data[i]->GetEfficiency(3, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
1643     effArray[i] = eff;
1644
1645     eff->GetXaxis()->SetRangeUser(0, 15);
1646     eff->GetYaxis()->SetRangeUser(0, 1.1);
1647     eff->SetStats(kFALSE);
1648     eff->SetTitle(";true multiplicity;Efficiency");
1649     eff->SetLineColor(colors[i]);
1650     eff->SetMarkerColor(colors[i]);
1651     eff->SetMarkerStyle(markers[i]);
1652
1653     if (i == 3)
1654     {
1655       for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1656         eff->SetBinError(bin, 0);
1657
1658       // loop over cross section combinations
1659       for (Int_t j=1; j<7; ++j)
1660       {
1661         AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
1662         mult->LoadHistograms(Form("Multiplicity_%d", j));
1663
1664         TH1* eff2 = mult->GetEfficiency(3, (AliMultiplicityCorrection::EventType) eventType);
1665
1666         for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1667         {
1668           // TODO we could also do assymetric errors here
1669           Float_t deviation = 10.0 * TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
1670
1671           eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), (Double_t) deviation));
1672         }
1673       }
1674
1675       for (Int_t bin=1; bin<=20; bin++)
1676         if (eff->GetBinContent(bin) > 0)
1677           Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1678
1679       effError = (TH1*) eff->Clone("effError");
1680       effError->Reset();
1681
1682       for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
1683         if (eff->GetBinContent(bin) > 0)
1684           effError->SetBinContent(bin, eff->GetBinError(bin) / eff->GetBinContent(bin));
1685
1686       effError->DrawCopy("SAME HIST");
1687     }
1688
1689     eff->SetBinContent(1, 0);
1690     eff->SetBinError(1, 0);
1691
1692     canvas->cd();
1693     if (i == 0)
1694     {
1695       eff->DrawCopy("P");
1696     }
1697     else
1698       eff->DrawCopy("SAME P");
1699
1700     legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
1701   }
1702
1703   legend->AddEntry(effError, "relative syst. uncertainty #times 10");
1704
1705   legend->Draw();
1706
1707   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1708 }
1709
1710 void ModelDependencyPlot()
1711 {
1712   gSystem->Load("libPWG0base");
1713
1714   TFile::Open("multiplicityMC_3M.root");
1715   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1716   mult->LoadHistograms("Multiplicity");
1717
1718   TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
1719
1720   TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
1721   canvas->SetGridx();
1722   canvas->SetGridy();
1723   //canvas->SetRightMargin(0.05);
1724   //canvas->SetTopMargin(0.05);
1725
1726   canvas->Divide(2, 1);
1727
1728   canvas->cd(2);
1729   gPad->SetLogy();
1730  
1731   Int_t selectedMult = 30;
1732   Int_t yMax = 200000;
1733
1734   TH1* full = proj->ProjectionX("full");
1735   TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult)); 
1736
1737   full->SetStats(kFALSE);
1738   full->GetXaxis()->SetRangeUser(0, 200);
1739   full->GetYaxis()->SetRangeUser(5, yMax);
1740   full->SetTitle(";multiplicity");
1741
1742   selected->SetLineColor(0);
1743   selected->SetMarkerColor(2);
1744   selected->SetMarkerStyle(7);
1745
1746   full->Draw();
1747   selected->Draw("SAME P");
1748
1749   TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
1750   legend->SetFillColor(0);
1751   legend->AddEntry(full, "true");
1752   legend->AddEntry(selected, "measured");
1753   legend->Draw();
1754  
1755   TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
1756   line->SetLineWidth(2);
1757   line->Draw();
1758
1759   canvas->cd(1);
1760   gPad->SetLogy();
1761
1762   full = proj->ProjectionY("full2");
1763   selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
1764
1765   full->SetStats(kFALSE);
1766   full->GetXaxis()->SetRangeUser(0, 200);
1767   full->GetYaxis()->SetRangeUser(5, yMax);
1768   full->SetTitle(";multiplicity");
1769
1770   full->SetLineColor(0);
1771   full->SetMarkerColor(2);
1772   full->SetMarkerStyle(7);
1773
1774   full->Draw("P");
1775   selected->Draw("SAME");
1776
1777   legend->Draw();
1778
1779   line = new TLine(selectedMult, 5, selectedMult, yMax);
1780   line->SetLineWidth(2);
1781   line->Draw();
1782
1783   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1784 }
1785
1786 void SystematicpTSpectrum()
1787 {
1788   gSystem->Load("libPWG0base");
1789
1790   TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
1791   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1792   mult->LoadHistograms("Multiplicity");
1793
1794   TFile::Open("multiplicityMC_100k_syst.root");
1795   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1796   mult2->LoadHistograms("Multiplicity");
1797
1798   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1799   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1800   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1801
1802   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1803   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1804
1805   DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
1806 }
1807
1808 // to be deleted
1809 /*void covMatrix(Bool_t mc = kTRUE)
1810 {
1811   gSystem->Load("libPWG0base");
1812
1813   TFile::Open(correctionFile);
1814   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1815   mult->LoadHistograms("Multiplicity");
1816
1817   TFile::Open(measuredFile);
1818   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1819   mult2->LoadHistograms("Multiplicity");
1820
1821   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1822
1823   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1824
1825   mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
1826 }*/
1827
1828 Double_t FitPtFunc(Double_t *x, Double_t *par)
1829 {
1830   Double_t xx = x[0];
1831
1832   Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
1833   Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
1834
1835   const Float_t kTransitionWidth = 0;
1836
1837   // power law part
1838   if (xx < par[0] - kTransitionWidth)
1839   {
1840     return val1;
1841   }
1842   /*else if (xx < par[0] + kTransitionWidth)
1843   {
1844     // smooth transition
1845     Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
1846     return (1 - factor) * val1 + factor * val2;
1847   }*/
1848   else
1849   {
1850     return val2;
1851   }
1852 }
1853
1854 void FitPt(const char* fileName = "firstplots100k_truept.root")
1855 {
1856   gSystem->Load("libPWG0base");
1857
1858   TFile::Open(fileName);
1859
1860   /*
1861   // merge corrections
1862   AliCorrection* correction[4];
1863   TList list;
1864
1865   for (Int_t i=0; i<4; ++i)
1866   {
1867     Printf("correction %d", i);
1868
1869     TString name; name.Form("correction_%d", i);
1870     correction[i] = new AliCorrection(name, name);
1871     correction[i]->LoadHistograms();
1872
1873     if (i > 0)
1874       list.Add(correction[i]);
1875   }
1876
1877   correction[0]->Merge(&list);
1878
1879   TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
1880
1881   // limit vtx, eta axis
1882   gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
1883   gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
1884
1885   TH1* genePt = gene->Project3D("z");*/
1886   TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
1887   genePt->Sumw2();
1888
1889   //genePt->Scale(1.0 / genePt->Integral());
1890
1891   // normalize by bin width
1892   for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1893     genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1894
1895   /// genePt->GetXaxis()->GetBinCenter(x));
1896
1897   genePt->GetXaxis()->SetRangeUser(0, 7.9);
1898   //genePt->GetYaxis()->SetTitle("a.u.");
1899
1900   TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
1901   //func->SetLineColor(2);
1902   func->SetParameters(1, -1, 1, 1, 1);
1903   func->SetParLimits(3, 1, 10);
1904   func->SetParLimits(4, 0, 10);
1905
1906   //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
1907
1908   //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
1909   //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
1910   //func->FixParameter(0, 0.314);
1911   //func->SetParLimits(0, 0.1, 0.3);
1912
1913   genePt->SetMarkerStyle(25);
1914   genePt->SetTitle("");
1915   genePt->SetStats(kFALSE);
1916   genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
1917   //func->Draw("SAME");
1918
1919   // fit only exp. part
1920   func->SetParameters(1, -1);
1921   func->FixParameter(2, 0);
1922   func->FixParameter(3, 1);
1923   func->FixParameter(4, 1);
1924   genePt->Fit(func, "0", "", 0.2, 1);
1925
1926   new TCanvas;
1927   genePt->DrawCopy("P");
1928   func->SetRange(0.02, 8);
1929   func->DrawCopy("SAME");
1930   gPad->SetLogy();
1931
1932   // now fix exp. parameters and fit second part
1933   Double_t param0 = func->GetParameter(0);
1934   Double_t param1 = func->GetParameter(1);
1935   func->SetParameters(0, -1, 1, 1, 1);
1936   func->FixParameter(0, 0);
1937   func->FixParameter(1, -1);
1938   func->ReleaseParameter(2);
1939   func->ReleaseParameter(3);
1940   func->ReleaseParameter(4);
1941   func->SetParLimits(3, 1, 10);
1942   func->SetParLimits(4, 0, 10);
1943
1944   genePt->Fit(func, "0", "", 1.5, 4);
1945
1946   new TCanvas;
1947   genePt->DrawCopy("P");
1948   func->SetRange(0.02, 8);
1949   func->DrawCopy("SAME");
1950   gPad->SetLogy();
1951
1952   // fit both
1953   func->SetParameter(0, param0);
1954   func->SetParameter(1, param1);
1955   func->ReleaseParameter(0);
1956   func->ReleaseParameter(1);
1957
1958   new TCanvas;
1959   genePt->DrawCopy("P");
1960   func->SetRange(0.02, 5);
1961   func->DrawCopy("SAME");
1962   gPad->SetLogy();
1963
1964   genePt->Fit(func, "0", "", 0.2, 4);
1965
1966   TCanvas* canvas = new TCanvas("FitPt", "FitPt", 800, 400);
1967   canvas->Divide(2, 1);
1968   canvas->cd(1);
1969
1970   gPad->SetGridx();
1971   gPad->SetGridy();
1972   gPad->SetLeftMargin(0.13);
1973   gPad->SetRightMargin(0.05);
1974   gPad->SetTopMargin(0.05);
1975
1976   genePt->GetXaxis()->SetRangeUser(0, 4.9);
1977   genePt->GetYaxis()->SetRangeUser(1e-2, 1e4);
1978   genePt->GetYaxis()->SetTitleOffset(1.4);
1979   genePt->GetXaxis()->SetTitleOffset(1.1);
1980   genePt->DrawCopy("P");
1981   func->SetRange(0.02, 5);
1982   func->DrawCopy("SAME");
1983   gPad->SetLogy();
1984
1985   canvas->cd(2);
1986
1987   TH1* genePtClone = (TH1*) genePt->Clone("genePtClone");
1988   genePtClone->Reset();
1989   genePtClone->DrawCopy("P");
1990
1991   gPad->SetGridx();
1992   gPad->SetGridy();
1993   gPad->SetLeftMargin(0.13);
1994   gPad->SetRightMargin(0.05);
1995   gPad->SetTopMargin(0.05);
1996
1997   func->DrawCopy("SAME");
1998   gPad->SetLogy();
1999
2000   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2001
2002   TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2003
2004   TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2005
2006   TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2007
2008   for (Int_t param=0; param<5; param++)
2009   {
2010     for (Int_t sign=0; sign<2; sign++)
2011     {
2012       TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign));  // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2013       func2->SetParameters(func->GetParameters());
2014       //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2015
2016       Float_t factor = ((sign == 0) ? 0.9 : 1.1);
2017       func2->SetParameter(param, func2->GetParameter(param) * factor);
2018       //func2->Print();
2019
2020       canvas->cd(2);
2021       func2->SetLineWidth(1);
2022       func2->SetLineColor(2);
2023       func2->DrawCopy("SAME");
2024
2025       canvas2->cd();
2026       TH1* second = func2->GetHistogram();
2027       second->Divide(first);
2028       second->SetLineColor(param + 1);
2029       second->GetYaxis()->SetRangeUser(0, 2);
2030       second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2031       second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2032     }
2033   }
2034
2035   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2036   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2037
2038   file->Close();
2039 }
2040
2041 void DrawSystematicpT()
2042 {
2043   TFile* file = TFile::Open("SystematicpT.root");
2044
2045   TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2046   TH1* result2 = (TH1*) file->Get("result_unity");
2047
2048   TH1* mcHist[12];
2049   TH1* result[12];
2050
2051   Int_t nParams = 5;
2052
2053   for (Int_t id=0; id<nParams*2; ++id)
2054   {
2055     mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2056     result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2057   }
2058
2059   DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2060
2061   //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2062
2063   DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2064
2065   //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2066
2067   // does not make sense: mc is different
2068   //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2069 }
2070
2071 void SystematicpT(Bool_t chi2 = 1)
2072 {
2073   gSystem->Load("libPWG0base");
2074
2075   TFile::Open("ptspectrum900.root");
2076   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2077   mult->LoadHistograms("Multiplicity");
2078
2079   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2080
2081   TH1* mcHist[12];
2082   TH1* result[12];
2083
2084   Int_t nParams = 5;
2085
2086   for (Int_t param=0; param<nParams; param++)
2087   {
2088     for (Int_t sign=0; sign<2; sign++)
2089     {
2090       // calculate result with systematic effect
2091       TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2092       mult2->LoadHistograms("Multiplicity");
2093
2094       mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2095
2096       if (chi2)
2097       {
2098         mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2099         mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2100       }
2101       else
2102         mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2103
2104       Int_t id = param * 2 + sign;
2105
2106       mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2107       result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2108
2109       TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2110       DrawResultRatio(mcHist[id], result[id], tmp);
2111     }
2112   }
2113
2114   // calculate normal result
2115   TFile::Open("ptspectrum100_1.root");
2116   mult2->LoadHistograms("Multiplicity");
2117   TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2118
2119   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2120
2121   if (chi2)
2122   {
2123     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2124   }
2125   else
2126     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2127
2128   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2129
2130   TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2131   mcHist2->Write();
2132   result2->Write();
2133   for (Int_t id=0; id<nParams*2; ++id)
2134   {
2135     mcHist[id]->Write();
2136     result[id]->Write();
2137   }
2138   file->Close();
2139
2140   DrawSystematicpT();
2141 }
2142
2143 void SystematicpTCutOff()
2144 {
2145   gSystem->Load("libPWG0base");
2146   SetTPC();
2147
2148   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2149   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2150   mult->LoadHistograms("Multiplicity");
2151
2152   TFile::Open(measuredFile);
2153   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2154   mult2->LoadHistograms("Multiplicity");
2155
2156   // "normal" result
2157   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2158   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2159   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2160
2161   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2162   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2163
2164   // change of pt spectrum (down)
2165   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2166   AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2167   mult3->LoadHistograms("Multiplicity");
2168
2169   mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2170   mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2171   mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2172
2173   TH1* result2 = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2174
2175   Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
2176 }