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