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