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   residualHist->Sumw2();
 933
 934   Float_t chi2 = 0;
 935   for (Int_t i=1; i<=displayRange+1; ++i)
 936   {
 937     if (measured->GetBinError(i) > 0)
 938     {
 939       residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
 940       residual->SetBinError(i, 1);
 941
 942       residualHist->Fill(residual->GetBinContent(i));
 943       chi2 += residual->GetBinContent(i) * residual->GetBinContent(i);
 944     }
 945     else
 946     {
 947       residual->SetBinContent(i, 0);
 948       residual->SetBinError(i, 0);
 949     }
 950   }
 951
 952   Printf("chi2 / ndf = %f / %d = %f", chi2, displayRange+1, chi2 / (displayRange+1));
 953
 954   residual->GetYaxis()->SetTitle("Residuals:   (1/e) (M - R  #otimes U)");
 955   residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
 956   residual->DrawCopy();
 957
 958   TLine* line = new TLine(-0.5, 0, displayRange + 0.5, 0);
 959   line->SetLineWidth(2);
 960   line->Draw();
 961
 962   pad3->cd();
 963   residualHist->SetStats(kFALSE);
 964   residualHist->SetLabelSize(0.08, "xy");
 965   residualHist->Fit("gaus");
 966   residualHist->Draw("HIST");
 967   residualHist->FindObject("gaus")->Draw("SAME");
 968
 969   canvas->Modified();
 970   canvas->SaveAs(canvas->GetName());
 971
 972   //const char* epsName2 = "proj.eps";
 973   //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
 974   //canvas->SetGridx();
 975   //canvas->SetGridy();
 976
 977   //canvas->SaveAs(canvas->GetName());
 978 }
 979
 980 void chi2FluctuationResult()
 981 {
 982   loadlibs();
 983
 984   TFile::Open("chi2_noregularization.root");
 985   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 986   mult->LoadHistograms("Multiplicity");
 987
 988   TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
 989
 990   mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
 991
 992   TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_1");
 993   ((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetRangeUser(0, displayRange);
 994   ((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->SetRangeUser(0, displayRange);
 995   //((TH1*) canvas->FindObject("proj"))->GetXaxis()->SetTitle(GetMultTitle());
 996   //((TH1*) canvas->FindObject("fCurrentESD"))->GetXaxis()->->SetTitle(GetMultTitle(etaRange, kFALSE));
 997   canvas->SaveAs("chi2FluctuationResult.eps");
 998 }
 999
1000 void DrawUnfolded(const char* fileName, const char* eps)
1001 {
1002   loadlibs();
1003
1004   TFile::Open(fileName);
1005
1006   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1007   mult->LoadHistograms("Multiplicity");
1008
1009   TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult->GetMultiplicityVtx(etaRange)->GetNbinsX());
1010   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1011
1012   DrawResultRatio(mcHist, result, eps);
1013 }
1014
1015 void minimizationInfluenceAlpha()
1016 {
1017   loadlibs();
1018
1019   TFile::Open(measuredFile);
1020   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1021   mult2->LoadHistograms("Multiplicity");
1022
1023   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, mult2->GetMultiplicityVtx(etaRange)->GetNbinsX());
1024   mcHist->Scale(1.0 / mcHist->Integral());
1025   mcHist->SetStats(kFALSE);
1026   mcHist->SetTitle(";True multiplicity n in |#eta| < 1.0;P(n)");
1027
1028   TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 350);
1029   canvas->Divide(3, 1, 0.005);
1030
1031   TFile::Open("chi2compare-influencealpha/EvaluateChi2Method.root");
1032
1033   TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_3.162278");
1034   TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_07_2_10000.000000");
1035   TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_13_2_10000000.000000");
1036
1037   /*mcHist->Rebin(2);  mcHist->Scale(0.5);
1038   hist1->Rebin(2);   hist1->Scale(0.5);
1039   hist2->Rebin(2);   hist2->Scale(0.5);
1040   hist3->Rebin(2);   hist3->Scale(0.5);*/
1041
1042   mcHist->GetXaxis()->SetRangeUser(0, displayRange);
1043   mcHist->SetLabelSize(0.06, "xy");
1044   mcHist->SetTitleSize(0.06, "xy");
1045   mcHist->GetYaxis()->SetTitleOffset(1.5);
1046
1047   canvas->cd(1);
1048
1049   gPad->SetLogy();
1050   gPad->SetRightMargin(0.03);
1051   gPad->SetLeftMargin(0.19);
1052   gPad->SetTopMargin(0.05);
1053   gPad->SetBottomMargin(0.13);
1054   gPad->SetGridx();
1055   gPad->SetGridy();
1056   mcHist->Draw();
1057   hist1->SetMarkerStyle(5);
1058   hist1->SetMarkerColor(2);
1059   hist1->Draw("SAME PE");
1060
1061   canvas->cd(2);
1062   gPad->SetRightMargin(0.03);
1063   gPad->SetLeftMargin(0.19);
1064   gPad->SetTopMargin(0.05);
1065   gPad->SetBottomMargin(0.13);
1066   gPad->SetLogy();
1067   gPad->SetGridx();
1068   gPad->SetGridy();
1069   mcHist->Draw();
1070   hist2->SetMarkerStyle(5);
1071   hist2->SetMarkerColor(2);
1072   hist2->Draw("SAME PE");
1073
1074   canvas->cd(3);
1075   gPad->SetRightMargin(0.03);
1076   gPad->SetLeftMargin(0.19);
1077   gPad->SetTopMargin(0.05);
1078   gPad->SetBottomMargin(0.13);
1079   gPad->SetLogy();
1080   gPad->SetGridx();
1081   gPad->SetGridy();
1082   mcHist->Draw();
1083   hist3->SetMarkerStyle(5);
1084   hist3->SetMarkerColor(2);
1085   hist3->Draw("SAME PE");
1086
1087   canvas->SaveAs("minimizationInfluenceAlpha.eps");
1088 }
1089
1090 void NBDFit()
1091 {
1092   gSystem->Load("libPWG0base");
1093
1094   TFile::Open(correctionFile);
1095   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1096   mult->LoadHistograms("Multiplicity");
1097
1098   TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
1099   fCurrentESD->Sumw2();
1100   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1101
1102   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])");
1103   func->SetParNames("scaling", "averagen", "k");
1104   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
1105   func->SetParLimits(1, 0.001, 1000);
1106   func->SetParLimits(2, 0.001, 1000);
1107   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
1108
1109   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
1110   lognormal->SetParNames("scaling", "mean", "sigma");
1111   lognormal->SetParameters(1, 1, 1);
1112   lognormal->SetParLimits(0, 0, 10);
1113   lognormal->SetParLimits(1, 0, 100);
1114   lognormal->SetParLimits(2, 1e-3, 10);
1115
1116   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
1117   fCurrentESD->SetStats(kFALSE);
1118   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
1119   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
1120   fCurrentESD->Draw("HIST");
1121   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1122   fCurrentESD->Fit(func, "W0", "", 0, 50);
1123   func->SetRange(0, 100);
1124   func->Draw("SAME");
1125   printf("chi2 = %f\n", func->GetChisquare());
1126
1127   fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
1128   lognormal->SetLineColor(2);
1129   lognormal->SetLineStyle(2);
1130   lognormal->SetRange(0, 100);
1131   lognormal->Draw("SAME");
1132
1133   canvas->SaveAs("NBDFit.eps");
1134 }
1135
1136 void StartingConditions()
1137 {
1138   // data generated by runMultiplicitySelector.C StartingConditions
1139
1140   const char* name = "StartingConditions";
1141
1142   TFile* file = TFile::Open(Form("%s.root", name));
1143
1144   TCanvas* canvas = new TCanvas(name, name, 1200, 600);
1145   canvas->Divide(2, 1);
1146
1147   TH1* mc = (TH1*) file->Get("mc");
1148   mc->Sumw2();
1149   mc->Scale(1.0 / mc->Integral(2, displayRange));
1150
1151   //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
1152
1153   TLegend* legend = new TLegend(0.6, 0.7, 0.99, 0.99);
1154   legend->SetFillColor(0);
1155   legend->SetTextSize(0.04);
1156
1157   const char* names[] = { "True", "Measured 1", "Measured 2", "Measured 3", "NBD", "Flat" };
1158
1159   for (Int_t i=0; i<6; ++i)
1160   {
1161     Int_t id = i;
1162     if (id > 2)
1163       id += 2;
1164
1165     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1166     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
1167
1168     chi2Result->Scale(1.0 / chi2Result->Integral(2, displayRange));
1169     bayesResult->Scale(1.0 / bayesResult->Integral(2, displayRange));
1170
1171     chi2Result->Divide(mc, chi2Result, 1, 1, "");
1172     bayesResult->Divide(mc, bayesResult, 1, 1, "");
1173
1174     chi2Result->SetTitle(Form("a) #chi^{2}-minimization;%s;MC / unfolded", GetMultLabel()));
1175     chi2Result->GetXaxis()->SetRangeUser(0, displayRange);
1176     chi2Result->GetYaxis()->SetRangeUser(0.7, 1.3);
1177     chi2Result->GetYaxis()->SetTitleOffset(1.7);
1178     //chi2Result->SetMarkerStyle(marker[i]);
1179     chi2Result->SetLineColor(i+1);
1180     chi2Result->SetMarkerColor(i+1);
1181     chi2Result->SetStats(kFALSE);
1182
1183     bayesResult->SetTitle(Form("b) Bayesian unfolding;%s;MC / unfolded", GetMultLabel()));
1184     bayesResult->GetXaxis()->SetRangeUser(0, displayRange);
1185     bayesResult->GetYaxis()->SetRangeUser(0.7, 1.3);
1186     bayesResult->GetYaxis()->SetTitleOffset(1.7);
1187     bayesResult->SetStats(kFALSE);
1188     //bayesResult->SetLineColor(2);
1189     bayesResult->SetLineColor(i+1);
1190
1191     canvas->cd(1);
1192     gPad->SetRightMargin(0.05);
1193     gPad->SetLeftMargin(0.12);
1194     gPad->SetGridx();
1195     gPad->SetGridy();
1196     chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1197
1198     canvas->cd(2);
1199     gPad->SetRightMargin(0.05);
1200     gPad->SetLeftMargin(0.12);
1201     gPad->SetGridx();
1202     gPad->SetGridy();
1203     bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1204
1205     //TLine* line = new TLine(0, 1, 150, 1);
1206     //line->Draw();
1207
1208     legend->AddEntry(chi2Result, names[i]);
1209   }
1210
1211   canvas->cd(1);
1212   legend->Draw();
1213   canvas->cd(2);
1214   legend->Draw();
1215   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1216 }
1217
1218 void StatisticsPlot()
1219 {
1220   const char* name = "StatisticsPlot";
1221
1222   TFile* file = TFile::Open(Form("%s.root", name));
1223
1224   TCanvas* canvas = new TCanvas(name, name, 600, 400);
1225
1226   TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
1227   fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1228   fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1229   fitResultsChi2->Draw("AP");
1230
1231   TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1232   fitResultsChi2->Fit(f);
1233
1234   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1235
1236   TH1* mc[5];
1237   TH1* result[5];
1238
1239   const char* plotname = "chi2Result";
1240
1241   name = "StatisticsPlotRatios";
1242   canvas = new TCanvas(name, name, 600, 400);
1243
1244   for (Int_t i=0; i<5; ++i)
1245   {
1246     mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1247     result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1248
1249     result[i]->SetLineColor(i+1);
1250     result[i]->Draw(((i == 0) ? "" : "SAME"));
1251   }
1252 }
1253
1254 void Draw2Unfolded(const char* file1, const char* file2, const char* output)
1255 {
1256   loadlibs();
1257
1258   TFile::Open(file1);
1259   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1260   mult->LoadHistograms("Multiplicity");
1261
1262   // result with systematic effect
1263   TFile::Open(file2);
1264   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1265   mult2->LoadHistograms("Multiplicity");
1266
1267   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1268   TH1* result1 = (TH1*) mult2->GetMultiplicityESDCorrected(etaRange)->Clone("result1"); // from file2 (with syst)
1269   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");  // from file1 (without syst)
1270
1271   DrawResultRatio(mcHist, result1, "tmp1.eps");
1272   DrawResultRatio(mcHist, result2, "tmp2.eps");
1273   Draw2ResultRatio(mcHist, result1, result2, output);
1274 }
1275
1276 void PythiaPhojet()
1277 {
1278   loadlibs();
1279
1280   displayRange = 55;
1281   Draw2Unfolded("self.root", "pythia.root", "test.eps");
1282
1283   canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last();
1284   pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First();
1285   pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last();
1286   legend = (TLegend*)pad1->GetListOfPrimitives()->Last();
1287
1288   ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (Pythia / Phojet)");
1289   ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (Pythia)");
1290   ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (Phojet)");
1291   canvas->SaveAs("PythiaPhojet.eps");
1292 }
1293
1294 void Misalignment()
1295 {
1296   loadlibs();
1297
1298   Draw2Unfolded("chi2_ideal.root", "chi2_misaligned.root", "test.eps");
1299
1300   canvas = (TCanvas*) gROOT->GetListOfCanvases()->Last();
1301   pad1 = (TCanvas*)canvas->GetListOfPrimitives()->First();
1302   pad2 = (TCanvas*)canvas->GetListOfPrimitives()->Last();
1303   legend = (TLegend*)pad1->GetListOfPrimitives()->Last();
1304
1305   ((TH1*)pad2->GetListOfPrimitives()->At(1))->GetYaxis()->SetTitle("Ratio (misaligned / realigned)");
1306   ((TLegendEntry*)legend->GetListOfPrimitives()->At(1))->SetLabel("Unfolded distribution (misaligned)");
1307   ((TLegendEntry*)legend->GetListOfPrimitives()->At(2))->SetLabel("Unfolded distribution (realigned)");
1308   canvas->SaveAs("SystematicMisalignment.eps");
1309 }
1310
1311 void SystematicLowEfficiency()
1312 {
1313   Draw2Unfolded("chi2.root", "chi2_loweff_5.root", "SystematicLowEfficiency.eps");
1314 }
1315
1316 void SystematicLowEfficiency2()
1317 {
1318   loadlibs();
1319
1320   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("chi2.root");
1321   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1322   TH1* mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1323   result2->Scale(1.0 / result2->Integral(2, displayRange));
1324   result2->Scale(mcHist->Integral(2, displayRange));
1325
1326   mult = AliMultiplicityCorrection::Open("chi2_loweff_5.root");
1327   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1328
1329   mult = AliMultiplicityCorrection::Open("chi2_loweff_2.root");
1330   TH1* ratio2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio2");
1331   ratio2->Scale(1.0 / ratio2->Integral(2, displayRange));
1332   ratio2->Scale(mcHist->Integral(2, displayRange));
1333   ratio2->Divide(result2);
1334
1335   mult = AliMultiplicityCorrection::Open("chi2_loweff_1.root");
1336   TH1* ratio3 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("ratio3");
1337   ratio3->Scale(1.0 / ratio3->Integral(2, displayRange));
1338   ratio3->Scale(mcHist->Integral(2, displayRange));
1339   ratio3->Divide(result2);
1340
1341   Draw2ResultRatios(mcHist, result1, result2, ratio2, ratio3, "SystematicLowEfficiency2.eps");
1342 }
1343
1344 void SystematicMisalignment()
1345 {
1346   gSystem->Load("libPWG0base");
1347
1348   TFile::Open(correctionFile);
1349   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1350   mult->LoadHistograms("Multiplicity");
1351
1352   TFile::Open("multiplicityMC_100k_fullmis.root");
1353   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1354   mult2->LoadHistograms("Multiplicity");
1355
1356   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1357   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1358   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1359
1360   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1361   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1362
1363   DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1364 }
1365
1366 void SystematicMisalignmentTPC()
1367 {
1368   gSystem->Load("libPWG0base");
1369
1370   SetTPC();
1371
1372   TFile::Open(correctionFile);
1373   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1374   mult->LoadHistograms("Multiplicity");
1375
1376   TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1377   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1378   mult2->LoadHistograms("Multiplicity");
1379
1380   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1381   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1382   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1383
1384   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1385   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1386
1387   DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1388 }
1389
1390 void LowMomentumEffectSPD()
1391 {
1392   // 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
1393   // only a normal acceptance region is considered to not get a bias by the edges
1394
1395   loadlibs();
1396   TFile::Open("multiplicity.root");
1397
1398
1399   AliCorrection* correction[8];
1400   Float_t values[3];
1401
1402   for (Int_t loop=0; loop<3; loop++)
1403   {
1404     Float_t sumGen = 0;
1405     Float_t sumMeas = 0;
1406
1407     Printf("loop %d", loop);
1408     for (Int_t i=0; i<4; ++i)
1409     {
1410       Printf("correction %d", i);
1411
1412       TString name; name.Form("correction_%d", i);
1413       correction[i] = new AliCorrection(name, name);
1414       correction[i]->LoadHistograms();
1415
1416       TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1417       TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1418
1419       Float_t vtxRange = 5.9;
1420       gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1421       meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1422
1423       Float_t etaRange = 0.99;
1424       gene->GetYaxis()->SetRangeUser(-etaRange, etaRange);
1425       meas->GetYaxis()->SetRangeUser(-etaRange, etaRange);
1426
1427       TH1* genePt = gene->Project3D(Form("z_%d", i));
1428       TH1* measPt = meas->Project3D(Form("z_%d", i));
1429
1430       if (loop > 0)
1431       {
1432         for (Int_t x=1; x<=genePt->GetNbinsX(); x++)
1433         {
1434           Float_t pt = genePt->GetXaxis()->GetBinCenter(x);
1435           //Printf("%f", pt);
1436           if (pt < 0.2)
1437           {
1438             Float_t factor = 1;
1439             if (loop == 1)
1440               factor = 1.5 - pt / 0.2 * 0.5;
1441             if (loop == 2)
1442               factor = 0.5 + pt / 0.2 * 0.5;
1443             //Printf("%f", factor);
1444             genePt->SetBinContent(x, genePt->GetBinContent(x) * factor);
1445             measPt->SetBinContent(x, measPt->GetBinContent(x) * factor);
1446           }
1447         }
1448       }
1449
1450       //new TCanvas; genePt->DrawCopy(); measPt->DrawCopy("SAME");
1451
1452       sumGen += genePt->Integral();
1453       sumMeas += measPt->Integral();
1454
1455       Float_t average = measPt->Integral() / genePt->Integral();
1456
1457       Printf("The average efficiency of this correction is %f", average);
1458     }
1459
1460     Float_t average = sumMeas / sumGen;
1461
1462     Printf("The average efficiency of all corrections is %f", average);
1463     values[loop] = average;
1464   }
1465
1466   Printf("relative is %f and %f", values[1] / values[0], values[2] / values[0]);
1467 }
1468
1469
1470 void EfficiencySpecies(Bool_t addDecayStopped = kFALSE)
1471 {
1472   loadlibs();
1473
1474   Int_t marker[] = {24, 25, 26, 27};
1475   Int_t color[] = {1, 2, 4, 3};
1476
1477   // SPD TPC
1478   //const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
1479   //const char* fileName[] = { "spd/multiplicity.root", "tpc/multiplicity.root" };
1480   const char* fileName[] = { "multiplicity.root", "multiplicity.root" };
1481   Float_t etaRangeArr[] = {0.49, 0.9};
1482   const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1483
1484   TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1485   canvas->Divide(2, 1);
1486
1487   TCanvas* canvas3 = new TCanvas("EfficiencySpecies_comb", "EfficiencySpecies_comb", 600, 600);
1488   gPad->SetGridx();
1489   gPad->SetGridy();
1490   gPad->SetRightMargin(0.05);
1491   gPad->SetTopMargin(0.05);
1492
1493   TLegend* legends[2];
1494
1495   for (Int_t loop=1; loop<2; ++loop)
1496   {
1497     Printf("%s", fileName[loop]);
1498
1499     TCanvas* canvas2 = new TCanvas(Form("EfficiencySpecies_%d", loop), Form("EfficiencySpecies_%d", loop), 600, 600);
1500     gPad->SetGridx();
1501     gPad->SetGridy();
1502     gPad->SetRightMargin(0.05);
1503     gPad->SetTopMargin(0.05);
1504
1505     AliCorrection* correction[8];
1506
1507     canvas->cd(loop+1);
1508
1509     gPad->SetGridx();
1510     gPad->SetGridy();
1511     gPad->SetRightMargin(0.05);
1512
1513     TLegend* legend = new TLegend(0.6, 0.4, 0.85, 0.6);
1514     legend->SetFillColor(0);
1515     legend->SetEntrySeparation(0.2);
1516     legend->SetTextSize(gStyle->GetTextSize());
1517
1518     legends[loop] = new TLegend(0.4+loop*0.3, 0.2, 0.6+loop*0.3, 0.5);
1519     legends[loop]->SetFillColor(0);
1520     legends[loop]->SetEntrySeparation(0.2);
1521     legends[loop]->SetTextSize(gStyle->GetTextSize());
1522     legends[loop]->SetHeader((loop == 0) ? "SPD" : "TPC");
1523
1524     Float_t below = 0;
1525     Float_t total = 0;
1526
1527     TFile* file = TFile::Open(fileName[loop]);
1528     if (!file)
1529     {
1530       Printf("Could not open %s", fileName[loop]);
1531       return;
1532     }
1533
1534     Float_t sumGen = 0;
1535     Float_t sumMeas = 0;
1536
1537     for (Int_t i=0; i<3; ++i)
1538     {
1539       Printf("correction %d", i);
1540
1541       TString name; name.Form("correction_%d", i);
1542       correction[i] = new AliCorrection(name, name);
1543       correction[i]->LoadHistograms();
1544
1545       TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1546       TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1547
1548       // limit vtx axis
1549       Float_t vtxRange = 3.9;
1550       gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1551       meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1552
1553       // 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)
1554       /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1555         for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1556         {
1557           gene->SetBinContent(x, 0, z, 0);
1558           gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1559           meas->SetBinContent(x, 0, z, 0);
1560           meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1561         }*/
1562
1563       // limit eta axis
1564       Float_t etaBegin = -etaRangeArr[loop];
1565       Float_t etaEnd   = etaRangeArr[loop];
1566       //etaBegin = 0.01;
1567       //etaEnd = -0.01;
1568       gene->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
1569       meas->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
1570
1571       TH1* genePt = gene->Project3D(Form("z_%d", i));
1572       TH1* measPt = meas->Project3D(Form("z_%d", i));
1573
1574       genePt->Sumw2();
1575       measPt->Sumw2();
1576
1577       for (Int_t x=0; x<=genePt->GetNbinsX()+1; x++)
1578       {
1579         genePt->SetBinError(x, TMath::Sqrt(genePt->GetBinContent(x)));
1580         measPt->SetBinError(x, TMath::Sqrt(measPt->GetBinContent(x)));
1581       }
1582
1583       sumGen += genePt->Integral();
1584       sumMeas += measPt->Integral();
1585
1586       TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1587       effPt->Reset();
1588       effPt->Divide(measPt, genePt, 1, 1, "B");
1589
1590       Int_t bin = 0;
1591       for (bin=20; bin>=1; bin--)
1592       {
1593         if (effPt->GetBinContent(bin) < 0.5)
1594           break;
1595       }
1596
1597       Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1598
1599       Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1600       Printf("%.4f of the particles are below that momentum", fraction);
1601
1602       below += genePt->Integral(1, bin);
1603       total += genePt->Integral();
1604
1605       effPt->SetLineColor(color[i]);
1606       effPt->SetMarkerColor(color[i]);
1607       effPt->SetMarkerStyle(marker[i]);
1608       effPt->SetMarkerSize(2);
1609
1610       effPt->GetXaxis()->SetRangeUser(0, 1);
1611       effPt->GetYaxis()->SetRangeUser(0.001, 1);
1612
1613       effPt->GetXaxis()->SetTitleOffset(1.1);
1614       effPt->GetYaxis()->SetTitleOffset(1.2);
1615
1616       effPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1617
1618       effPt->SetStats(kFALSE);
1619       effPt->SetTitle(titles[loop]);
1620       effPt->GetYaxis()->SetTitle("Efficiency");
1621
1622       canvas->cd(loop+1);
1623       effPt->DrawCopy((i == 0) ? "" : "SAME");
1624
1625       canvas2->cd();
1626       effPt->SetTitle("");
1627       effPt->DrawCopy((i == 0) ? "" : "SAME");
1628
1629       canvas3->cd();
1630       effPtClone = (TH1*) effPt->Clone("effPtClone");
1631       effPtClone->SetMarkerStyle(marker[i]-4*loop);
1632       effPtClone->DrawCopy((i == 0 && loop == 0) ? "" : "SAME");
1633
1634       legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1635       legends[loop]->AddEntry(effPtClone, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1636       //legend2->AddEntry(effPt, Form("%s %s", (loop == 0) ? "SPD" : "TPC", ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}"))), "P");
1637
1638       if (addDecayStopped)
1639       {
1640         name.Form("correction_%d", i+4);
1641         corr = new AliCorrection(name, name);
1642         corr->LoadHistograms();
1643
1644         TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
1645         TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
1646
1647         // limit axes
1648         gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1649         meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1650         gene->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1651         meas->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1652
1653         TH1* decayed = gene->Project3D(Form("z_%d", i+4));
1654         TH1* stopped = meas->Project3D(Form("z_%d", i+4));
1655
1656         Printf("%d: %d decayed, %d stopped, out of %d", i, (Int_t) decayed->Integral(), (Int_t) stopped->Integral(), (Int_t) genePt->Integral());
1657
1658         decayed->Divide(decayed, genePt, 1, 1, "B");
1659         stopped->Divide(stopped, genePt, 1, 1, "B");
1660
1661         decayed->SetMarkerStyle(20);
1662         stopped->SetMarkerStyle(21);
1663         stopped->SetMarkerColor(2);
1664
1665         new TCanvas(Form("all_%d_%d", loop, i), Form("all_%d_%d", loop, i), 600, 600);
1666         effPt->DrawCopy();
1667         decayed->DrawCopy("SAME");
1668         stopped->DrawCopy("SAME");
1669
1670         decayed->Add(stopped);
1671         decayed->Add(effPt);
1672         decayed->SetMarkerStyle(22);
1673         decayed->SetMarkerColor(4);
1674         decayed->DrawCopy("SAME");
1675       }
1676
1677     }
1678
1679     Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1680
1681     Printf("%f measured, %f generated, effiency: %f", sumGen, sumMeas, sumMeas / sumGen);
1682
1683     canvas->cd(loop+1);
1684     legend->Draw();
1685
1686     canvas2->cd();
1687     legend->Draw();
1688     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
1689   }
1690
1691   canvas->cd();
1692   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1693
1694   canvas3->cd();
1695   legends[0]->Draw();
1696   legends[1]->Draw();
1697   canvas3->SaveAs(Form("%s.eps", canvas3->GetName()));
1698 }
1699
1700 void DrawpTCutOff()
1701 {
1702 /*
1703 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt_ref.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1704 mv unfolded.root chi2_ptref.root
1705 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1706 mv unfolded.root chi2_pt0.root
1707 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1708 mv unfolded.root chi2_pt1.root
1709 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1710 mv unfolded.root chi2_pt0_25.root
1711 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1712 mv unfolded.root chi2_pt1_25.root
1713 */
1714
1715   loadlibs();
1716
1717   TH1* results[10];
1718   const char* files[] = { "chi2_ptref.root", "chi2_pt0.root", "chi2_pt1.root", "chi2_pt0_25.root", "chi2_pt1_25.root"};
1719
1720   Int_t nMax = 5;
1721   for (Int_t i = 0; i<nMax; ++i)
1722   {
1723     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(files[i]);
1724     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1725   }
1726
1727   const char* legendStrings[] = { "Reduced 50%", "Enhanced 50%", "Reduced 25%", "Enhanced 25%" };
1728   DrawRatio(results[0], nMax-1, results+1, "LowMomentumSyst.eps", kFALSE, legendStrings);
1729 }
1730
1731 void ParticleSpeciesComparison()
1732 {
1733   loadlibs();
1734
1735   TH1* results[10];
1736   TH1* mc = 0;
1737
1738   // loop over cases (normal, enhanced/reduced ratios)
1739   Int_t nMax = 9;
1740   for (Int_t i = 0; i<nMax; ++i)
1741   {
1742     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(Form("chi2_species_%d.root", i), Form("Multiplicity_%d", i));
1743     if (i == 0)
1744       mc = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymchist", 1, 1);
1745     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1746   }
1747
1748   DrawResultRatio(mc, results[0], "ParticleSpeciesComparison1_1.eps");
1749
1750   for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1751   {
1752     results[0]->SetBinError(i, 0);
1753     mc->SetBinError(i, 0);
1754   }
1755
1756   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" };
1757
1758   DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1759
1760   //not valid: draw chi2 uncertainty on top!
1761   /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1762   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1763   errorHist->SetLineColor(1);
1764   errorHist->SetLineWidth(2);
1765   TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1766   for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1767   {
1768     errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1769     errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1770   }
1771
1772   errorHist->DrawCopy("SAME");
1773   errorHist2->DrawCopy("SAME");*/
1774
1775   //canvas->SaveAs(canvas->GetName());
1776
1777   DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1778
1779   //errorHist->DrawCopy("SAME");
1780   //errorHist2->DrawCopy("SAME");
1781
1782   //canvas2->SaveAs(canvas2->GetName());
1783 }
1784
1785 /*void ParticleSpeciesComparison2()
1786 {
1787   gSystem->Load("libPWG0base");
1788
1789   const char* fileNameMC = "multiplicityMC_400k_syst.root";
1790   const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1791   Bool_t chi2 = 0;
1792
1793   TFile::Open(fileNameMC);
1794   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1795   mult->LoadHistograms();
1796
1797   TH1* mc[10];
1798   TH1* results[10];
1799
1800   // loop over cases (normal, enhanced/reduced ratios)
1801   Int_t nMax = 7;
1802   for (Int_t i = 0; i<nMax; ++i)
1803   {
1804     TString folder;
1805     folder.Form("Multiplicity_%d", i);
1806
1807     AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1808
1809     TFile::Open(fileNameESD);
1810     mult2->LoadHistograms();
1811
1812     mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1813
1814     if (chi2)
1815     {
1816       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1817       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1818       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1819     }
1820     else
1821     {
1822       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1823       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1824     }
1825
1826     //Float_t averageRatio = 0;
1827     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1828
1829     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1830
1831     TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1832     mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1833
1834     //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1835     //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1836
1837     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1838   }
1839
1840   DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1841 }*/
1842
1843 TH1* Invert(TH1* eff)
1844 {
1845   // calculate corr = 1 / eff
1846
1847   TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1848   corr->Reset();
1849
1850   for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1851   {
1852     if (eff->GetBinContent(i) > 0)
1853     {
1854       corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1855       corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1856     }
1857   }
1858
1859   return corr;
1860 }
1861
1862 void TriggerVertexCorrection()
1863 {
1864   //
1865   // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1866   //
1867
1868   loadlibs();
1869
1870   TFile::Open(correctionFile);
1871   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1872   mult->LoadHistograms("Multiplicity");
1873
1874   TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1875   TH1* corrNSD = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kNSD));
1876   TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1877
1878   TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 500);
1879   gPad->SetGridx();
1880   gPad->SetGridy();
1881   gPad->SetTopMargin(0.05);
1882   gPad->SetRightMargin(0.05);
1883
1884   corrINEL->SetStats(kFALSE);
1885   corrINEL->GetXaxis()->SetRangeUser(0, 12);
1886   corrINEL->GetYaxis()->SetRangeUser(0, 8);
1887   corrINEL->SetTitle(Form(";%s;Correction factor", GetMultLabel()));
1888   corrINEL->SetMarkerStyle(22);
1889   corrINEL->Draw("PE");
1890
1891   corrMB->SetStats(kFALSE);
1892   corrMB->SetLineColor(2);
1893   corrMB->SetMarkerStyle(25);
1894   corrMB->SetMarkerColor(2);
1895   corrMB->Draw("SAME PE");
1896
1897   corrNSD->SetLineColor(4);
1898   corrNSD->SetMarkerStyle(24);
1899   corrNSD->SetMarkerColor(4);
1900   corrNSD->Draw("SAME PE");
1901
1902   Printf("       MB  INEL  NSD");
1903   Printf("bin 0: %f %f %f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1));
1904   Printf("bin 1: %f %f %f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2));
1905
1906   TLegend* legend = new TLegend(0.3, 0.6, 0.85, 0.85);
1907   legend->SetFillColor(0);
1908   legend->AddEntry(corrINEL, "Correction to inelastic sample");
1909   legend->AddEntry(corrNSD, "Correction to NSD sample");
1910   legend->AddEntry(corrMB, "Correction to triggered sample");
1911   legend->SetTextSize(0.04);
1912
1913   legend->Draw();
1914
1915   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1916 }
1917
1918 void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1919 {
1920   loadlibs();
1921
1922   TFile::Open(correctionFile);
1923   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1924   mult->LoadHistograms("Multiplicity");
1925
1926   TFile::Open(measuredFile);
1927   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1928   mult2->LoadHistograms("Multiplicity");
1929
1930   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1931
1932   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1933
1934   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 1e5);
1935
1936   TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1937
1938   return;
1939
1940   TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1941
1942   TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1943
1944   if (!mc)
1945   {
1946     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1947     DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1948   }
1949
1950   TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
1951   errorResponse->Write();
1952   errorMeasured->Write();
1953   errorBoth->Write();
1954   file->Close();
1955 }
1956
1957 void DrawStatisticalUncertainty()
1958 {
1959   TFile::Open("StatisticalUncertainty.root");
1960
1961   errorResponse = (TH1*) gFile->Get("errorResponse");
1962   errorMeasured = (TH1*) gFile->Get("errorMeasured");
1963   errorBoth = (TH1*) gFile->Get("errorBoth");
1964
1965   TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1966   canvas->SetGridx();
1967   canvas->SetGridy();
1968   canvas->SetRightMargin(0.05);
1969   canvas->SetTopMargin(0.05);
1970
1971   errorResponse->SetLineColor(1);
1972   errorResponse->GetXaxis()->SetRangeUser(0, longDisplayRange);
1973   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1974   errorResponse->SetStats(kFALSE);
1975   errorResponse->SetTitle(";true multiplicity;Uncertainty");
1976
1977   errorResponse->Draw();
1978
1979   errorMeasured->SetLineColor(2);
1980   errorMeasured->Draw("SAME");
1981
1982   errorBoth->SetLineColor(4);
1983   errorBoth->Draw("SAME");
1984
1985   Printf("Average errorResponse: %f", errorResponse->Integral(2, displayRange) / (displayRange - 1));
1986   Printf("Average errorMeasured: %f", errorMeasured->Integral(2, displayRange) /  (displayRange - 1));
1987   Printf("Average errorBoth: %f", errorBoth->Integral(2, displayRange) /  (displayRange - 1));
1988
1989   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1990 }
1991
1992 void StatisticalUncertaintyCompare(const char* det = "SPD")
1993 {
1994   TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1995   TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1996   TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1997   TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1998
1999   TString str;
2000   str.Form("StatisticalUncertaintyCompare%s", det);
2001
2002   TCanvas* canvas = new TCanvas(str, str, 800, 500);
2003   canvas->SetGridx();
2004   canvas->SetGridy();
2005   canvas->SetRightMargin(0.05);
2006   canvas->SetTopMargin(0.05);
2007
2008   errorResponse->Scale(1.0 / sqrt(2));
2009   errorMeasured->Scale(1.0 / sqrt(2));
2010   errorBoth->Scale(1.0 / sqrt(2));
2011
2012   errorResponse->SetLineColor(1);
2013   errorResponse->GetXaxis()->SetRangeUser(0, displayRange);
2014   errorResponse->GetYaxis()->SetRangeUser(0, 0.18);
2015   errorResponse->SetStats(kFALSE);
2016   errorResponse->GetYaxis()->SetTitleOffset(1.2);
2017   errorResponse->SetTitle(Form(";%s;#sqrt{2}^{-1} #sigma(unfolded - unfolded_{0}) / unfolded_{0}", GetMultLabel()));
2018
2019   errorResponse->Draw();
2020
2021   errorMeasured->SetLineColor(2);
2022   errorMeasured->Draw("SAME");
2023
2024   errorBoth->SetLineColor(4);
2025   errorBoth->Draw("SAME");
2026
2027   TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
2028   TH1* errorResponse2 = (TH1*) file2->Get("errorResponse");
2029   TH1* errorMeasured2 = (TH1*) file2->Get("errorMeasured");
2030   TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
2031
2032   errorResponse2->Scale(1.0 / sqrt(2));
2033   errorMeasured2->Scale(1.0 / sqrt(2));
2034   errorBoth2->Scale(1.0 / sqrt(2));
2035
2036   errorResponse2->SetLineStyle(2);
2037   errorResponse2->Draw("SAME");
2038
2039   errorMeasured2->SetLineColor(2);
2040   errorMeasured2->SetLineStyle(2);
2041   errorMeasured2->Draw("SAME");
2042
2043   errorBoth2->SetLineColor(4);
2044   errorBoth2->SetLineStyle(2);
2045   errorBoth2->Draw("SAME");
2046
2047   TLegend* legend = new TLegend(0.2, 0.5, 0.8, 0.9);
2048   legend->SetFillColor(0);
2049   legend->SetTextSize(0.04);
2050   legend->AddEntry(errorBoth, "Both (Bayesian unfolding)");
2051   legend->AddEntry(errorMeasured, "Measured (Bayesian unfolding)");
2052   legend->AddEntry(errorResponse, "Response matrix (Bayesian unfolding)");
2053   legend->AddEntry(errorBoth2, "Both (#chi^{2}-minimization)");
2054   legend->AddEntry(errorMeasured2, "Measured (#chi^{2}-minimization)");
2055   legend->AddEntry(errorResponse2, "Response matrix (#chi^{2}-minimization)");
2056   legend->Draw();
2057
2058   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2059 }
2060
2061 void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
2062 {
2063  const char* files[] = { "multiplicityMC_nd.root", "multiplicityMC_sd.root", "multiplicityMC_dd.root", "multiplicityMC_xsection.root" };
2064
2065   loadlibs();
2066
2067   TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 600);
2068   canvas->SetGridx();
2069   canvas->SetGridy();
2070   canvas->SetRightMargin(0.05);
2071   canvas->SetTopMargin(0.05);
2072
2073   AliMultiplicityCorrection* data[4];
2074   TH1* effArray[4];
2075   TH1* effErrorArray[2];
2076
2077   Int_t markers[] = { 24, 25, 26, 5 };
2078   //Int_t markers[] = { 2, 25, 24, 5 };
2079   Int_t colors[] = { 1, 2, 4, 6 };
2080   //Int_t colors[] = { 1, 1, 1, 1 };
2081
2082   //TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
2083   TLegend* legend = new TLegend(0.3, 0.3, 0.9, 0.6);
2084   legend->SetTextSize(0.04);
2085   legend->SetFillColor(0);
2086
2087   for (Int_t i=0; i<4; ++i)
2088   {
2089     TString name;
2090     name.Form("Multiplicity_%d", i);
2091
2092     TFile::Open(files[i]);
2093     data[i] = new AliMultiplicityCorrection(name, name);
2094
2095     if (i < 3)
2096     {
2097       data[i]->LoadHistograms("Multiplicity");
2098     }
2099     else
2100       data[i]->LoadHistograms("Multiplicity_0");
2101
2102     TH1* eff = 0;
2103     if (eventType == -1)
2104     {
2105       eff = (TH1*) data[i]->GetTriggerEfficiency(etaRange)->Clone(Form("eff_%d", i));
2106     }
2107     else
2108       eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
2109     effArray[i] = eff;
2110
2111     eff->GetXaxis()->SetRangeUser(0, 15);
2112     eff->GetYaxis()->SetRangeUser(0, 1.19);
2113     eff->SetStats(kFALSE);
2114     eff->GetXaxis()->SetTitle(GetMultLabel());
2115     eff->GetYaxis()->SetTitle("Efficiency");
2116     eff->SetTitle("");
2117     eff->SetLineColor(colors[i]);
2118     eff->SetMarkerColor(colors[i]);
2119     eff->SetMarkerStyle(markers[i]);
2120
2121     if (i == 3)
2122     {
2123       // once for INEL, once for NSD
2124       for (AliMultiplicityCorrection::EventType eventType2 = AliMultiplicityCorrection::kINEL; eventType2 <= AliMultiplicityCorrection::kNSD; eventType2++)
2125       {
2126         effDiff = (TH1*) data[i]->GetEfficiency(etaRange, eventType2)->Clone(Form("effDiff_%d", i));
2127
2128         for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2129           effDiff->SetBinError(bin, 0);
2130
2131         // loop over cross section combinations
2132         for (Int_t j=1; j<7; ++j)
2133         {
2134           AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
2135           mult->LoadHistograms(Form("Multiplicity_%d", j));
2136
2137           TH1* eff2 = mult->GetEfficiency(etaRange, eventType2);
2138
2139           for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2140           {
2141             // TODO we could also do asymmetric errors here
2142             Float_t deviation = TMath::Abs(effDiff->GetBinContent(bin) - eff2->GetBinContent(bin));
2143
2144             effDiff->SetBinError(bin, TMath::Max(effDiff->GetBinError(bin), (Double_t) deviation));
2145           }
2146         }
2147
2148         for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2149         {
2150           //if (eventType2 == AliMultiplicityCorrection::kINEL)
2151             //eff->SetBinError(bin, 0);
2152             //eff->SetBinError(bin, effDiff->GetBinError(bin));
2153           if (bin < 20 && effDiff->GetBinContent(bin) > 0)
2154             Printf("Bin %d: Error: %.2f", bin, 100.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2155         }
2156
2157         if (uncertainty) {
2158                 TH1* effError = (TH1*) effDiff->Clone(Form("effError_%s", (eventType2 == AliMultiplicityCorrection::kINEL) ? "INEL" : "NSD"));
2159           effErrorArray[eventType2 - AliMultiplicityCorrection::kINEL] = effError;
2160                 effError->Reset();
2161
2162                 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2163                   if (effDiff->GetBinContent(bin) > 0)
2164                     effError->SetBinContent(bin, 1.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2165
2166                 effError->SetLineColor(1);
2167                 effError->SetMarkerStyle(1);
2168
2169           if (eventType2 == AliMultiplicityCorrection::kNSD)
2170             effError->SetLineStyle(2);
2171
2172           effError->DrawCopy("SAME HIST");
2173         }
2174       }
2175     }
2176
2177     canvas->cd();
2178     if (i == 0)
2179     {
2180       eff->DrawCopy("P");
2181     }
2182     else
2183       eff->DrawCopy("SAME P");
2184
2185     legend->AddEntry(eff, (((i == 0) ? "Non-diffractive" : ((i == 1) ? "Single-diffractive" : ((i == 2) ? "Double-diffractive" : "Pythia combined")))));
2186   }
2187
2188   if (uncertainty)
2189   {
2190     legend->AddEntry(effErrorArray[0], "Relative syst. uncertainty: inelastic");
2191     legend->AddEntry(effErrorArray[1], "Relative syst. uncertainty: NSD");
2192
2193     file = TFile::Open("uncertainty_xsection.root", "RECREATE");
2194     effErrorArray[0]->Write();
2195     effErrorArray[1]->Write();
2196     file->Close();
2197   }
2198
2199   legend->Draw();
2200
2201   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2202 }
2203
2204 void ModelDependencyPlot()
2205 {
2206   loadlibs();
2207
2208   TFile::Open("multiplicityMC.root");
2209   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2210   mult->LoadHistograms("Multiplicity");
2211
2212   hist = mult->GetCorrelation(etaRange);
2213
2214   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
2215   {
2216     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
2217     {
2218       hist->SetBinContent(0, y, z, 0);
2219       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2220     }
2221   }
2222
2223   TH2* proj = (TH2*) hist->Project3D("zy");
2224
2225   TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 1200, 600);
2226
2227   canvas->Divide(2, 1);
2228
2229   canvas->cd(2);
2230   gPad->SetLogy();
2231   gPad->SetGridx();
2232   gPad->SetGridy();
2233   gPad->SetTopMargin(0.05);
2234   gPad->SetRightMargin(0.05);
2235
2236   Int_t selectedMult = 30;
2237   Int_t yMax = 9e4;
2238
2239   TH1* full = proj->ProjectionX("full");
2240   TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
2241
2242   full->SetStats(kFALSE);
2243   full->GetXaxis()->SetRangeUser(0, displayRange);
2244   full->GetYaxis()->SetRangeUser(5, yMax);
2245   full->SetTitle(";Multiplicity;Entries");
2246
2247   selected->SetLineColor(0);
2248   selected->SetMarkerColor(2);
2249   selected->SetMarkerStyle(5);
2250
2251   full->Draw();
2252   selected->Draw("SAME P");
2253
2254   TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
2255   legend->SetFillColor(0);
2256   legend->AddEntry(full, "True");
2257   legend->AddEntry(selected, "Measured");
2258   legend->Draw();
2259
2260   TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
2261   line->SetLineWidth(2);
2262   line->Draw();
2263
2264   canvas->cd(1);
2265   gPad->SetLogy();
2266   gPad->SetGridx();
2267   gPad->SetGridy();
2268   gPad->SetTopMargin(0.05);
2269   gPad->SetRightMargin(0.05);
2270
2271   full = proj->ProjectionY("full2");
2272   selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
2273
2274   full->SetStats(kFALSE);
2275   full->GetXaxis()->SetRangeUser(0, displayRange);
2276   full->GetYaxis()->SetRangeUser(5, yMax);
2277   full->SetTitle(";Multiplicity;Entries");
2278
2279   full->SetLineColor(0);
2280   full->SetMarkerColor(2);
2281   full->SetMarkerStyle(5);
2282
2283   full->Draw("P");
2284   selected->Draw("SAME");
2285
2286   legend->Draw();
2287
2288   line = new TLine(selectedMult, 5, selectedMult, yMax);
2289   line->SetLineWidth(2);
2290   line->Draw();
2291
2292   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2293 }
2294
2295 void SystematicpTSpectrum()
2296 {
2297   gSystem->Load("libPWG0base");
2298
2299   TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
2300   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2301   mult->LoadHistograms("Multiplicity");
2302
2303   TFile::Open("multiplicityMC_100k_syst.root");
2304   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2305   mult2->LoadHistograms("Multiplicity");
2306
2307   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2308   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2309   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2310
2311   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2312   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2313
2314   DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
2315 }
2316
2317 void FitPt(const char* fileName)
2318 {
2319   // needs a MC file from the dNdEta analysis
2320
2321   TFile::Open(fileName);
2322
2323   TH1* genePt = (TH1*) gFile->Get("fdNdpT");
2324
2325   genePt->SetTitle(";p_{T} (GeV/c);1/p_{T} dN_{ch}/dp_{T} (GeV/c)^{-2}");
2326   // number of events
2327   genePt->Scale(1.0 / 287800);
2328   // bin width
2329   genePt->Scale(1.0 / genePt->GetXaxis()->GetBinWidth(1));
2330
2331   genePt->GetXaxis()->SetRangeUser(0, 0.4);
2332
2333   TF1* func = new TF1("func", "[1]*x*exp(x*[0])");
2334   func->SetParameters(-1, 1);
2335
2336   genePt->SetMarkerStyle(25);
2337   genePt->SetTitle("");
2338   genePt->SetStats(kFALSE);
2339
2340   new TCanvas;
2341   genePt->DrawCopy("P");
2342   func->DrawCopy("SAME");
2343   gPad->SetLogy();
2344
2345   genePt->Fit(func, "0", "", 0, 0.25);
2346   genePt->Fit(func, "0", "", 0, 0.25);
2347
2348   TCanvas* canvas = new TCanvas("FitPt", "FitPt", 600, 600);
2349
2350   gPad->SetGridx();
2351   gPad->SetGridy();
2352   gPad->SetLeftMargin(0.13);
2353   gPad->SetRightMargin(0.05);
2354   gPad->SetTopMargin(0.05);
2355
2356   //genePt->GetXaxis()->SetRangeUser(0, 1);
2357   genePt->GetYaxis()->SetRangeUser(2, 200);
2358   genePt->GetYaxis()->SetTitleOffset(1.4);
2359   genePt->GetXaxis()->SetTitleOffset(1.1);
2360   genePt->DrawCopy("P");
2361   //func->SetRange(0, 0.3);
2362   func->DrawCopy("SAME");
2363   gPad->SetLogy();
2364
2365   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2366
2367   TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2368
2369   TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2370
2371   TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2372
2373   for (Int_t param=0; param<2; param++)
2374   {
2375     for (Int_t sign=0; sign<2; sign++)
2376     {
2377       TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign));  // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2378       func2->SetParameters(func->GetParameters());
2379       //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2380
2381       Float_t factor = ((sign == 0) ? 0.75 : 1.25);
2382       func2->SetParameter(param, func2->GetParameter(param) * factor);
2383       //func2->Print();
2384
2385       canvas->cd();
2386       func2->SetLineWidth(2);
2387       func2->SetLineColor(2);
2388       func2->SetLineStyle(2);
2389       func2->DrawCopy("SAME");
2390
2391       canvas2->cd();
2392       TH1* second = func2->GetHistogram();
2393       second->Divide(first);
2394       second->SetLineColor(param + 1);
2395       // set to 1 above 0.2 GeV
2396       for (Int_t bin=second->GetXaxis()->FindBin(0.20001); bin<=second->GetNbinsX(); bin++)
2397         second->SetBinContent(bin, 1);
2398       second->GetYaxis()->SetRangeUser(0, 2);
2399       second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2400       second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2401     }
2402   }
2403
2404   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2405   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2406
2407   file->Close();
2408 }
2409
2410 void DrawSystematicpT()
2411 {
2412   TFile* file = TFile::Open("SystematicpT.root");
2413
2414   TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2415   TH1* result2 = (TH1*) file->Get("result_unity");
2416
2417   TH1* mcHist[12];
2418   TH1* result[12];
2419
2420   Int_t nParams = 5;
2421
2422   for (Int_t id=0; id<nParams*2; ++id)
2423   {
2424     mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2425     result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2426   }
2427
2428   DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2429
2430   //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2431
2432   DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2433
2434   //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2435
2436   // does not make sense: mc is different
2437   //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2438 }
2439
2440 void SystematicpT(Bool_t chi2 = 1)
2441 {
2442   gSystem->Load("libPWG0base");
2443
2444   TFile::Open("ptspectrum900.root");
2445   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2446   mult->LoadHistograms("Multiplicity");
2447
2448   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2449
2450   TH1* mcHist[12];
2451   TH1* result[12];
2452
2453   Int_t nParams = 5;
2454
2455   for (Int_t param=0; param<nParams; param++)
2456   {
2457     for (Int_t sign=0; sign<2; sign++)
2458     {
2459       // calculate result with systematic effect
2460       TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2461       mult2->LoadHistograms("Multiplicity");
2462
2463       mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2464
2465       if (chi2)
2466       {
2467         mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2468         mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2469       }
2470       else
2471         mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2472
2473       Int_t id = param * 2 + sign;
2474
2475       mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2476       result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2477
2478       TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2479       DrawResultRatio(mcHist[id], result[id], tmp);
2480     }
2481   }
2482
2483   // calculate normal result
2484   TFile::Open("ptspectrum100_1.root");
2485   mult2->LoadHistograms("Multiplicity");
2486   TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2487
2488   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2489
2490   if (chi2)
2491   {
2492     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2493   }
2494   else
2495     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2496
2497   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2498
2499   TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2500   mcHist2->Write();
2501   result2->Write();
2502   for (Int_t id=0; id<nParams*2; ++id)
2503   {
2504     mcHist[id]->Write();
2505     result[id]->Write();
2506   }
2507   file->Close();
2508
2509   DrawSystematicpT();
2510 }
2511
2512 void DrawSystematicpT2()
2513 {
2514   //displayRange = 200;
2515
2516   // read from file
2517   TFile* file = TFile::Open("SystematicpT2.root");
2518   TH1* mcHist = (TH1*) file->Get("mymc");
2519   TH1* result[12];
2520   result[0] = (TH1*) file->Get("result_unity");
2521   Int_t nParams = 5;
2522   for (Int_t id=0; id<nParams*2; ++id)
2523     result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2524
2525   DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2526   DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2527   DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2528 }
2529
2530 void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2531 {
2532   gSystem->Load("libPWG0base");
2533
2534   if (tpc)
2535   {
2536     SetTPC();
2537     TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2538   }
2539   else
2540     TFile::Open("ptspectrum100_1.root");
2541
2542   AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2543   measured->LoadHistograms("Multiplicity");
2544   TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2545
2546   TH1* result[12];
2547
2548   Int_t nParams = 5;
2549
2550   // -1 = unity change, 0...4 parameters
2551   for (Int_t id=-1; id<nParams*2; id++)
2552   {
2553     Int_t param = id / 2;
2554     Int_t sign = id % 2;
2555
2556     TString idStr;
2557     if (id == -1)
2558     {
2559       idStr = "unity";
2560     }
2561     else
2562       idStr.Form("%d_%d", param, sign);
2563
2564     // calculate result with systematic effect
2565     if (tpc)
2566     {
2567       TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2568     }
2569     else
2570       TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2571
2572     AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2573     response->LoadHistograms("Multiplicity");
2574
2575     response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2576
2577     if (chi2)
2578     {
2579       response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2580       response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2581     }
2582     else
2583       response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2584
2585     result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2586
2587     TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2588     DrawResultRatio(mcHist, result[id+1], tmp);
2589   }
2590
2591   TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2592   mcHist->Write();
2593   for (Int_t id=0; id<nParams*2+1; ++id)
2594     result[id]->Write();
2595   file->Close();
2596
2597   DrawSystematicpT2();
2598 }
2599
2600 void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2601 {
2602   // only needed for TPC
2603   SetTPC();
2604
2605   gSystem->Load("libPWG0base");
2606
2607   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2608   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2609   mult->LoadHistograms("Multiplicity");
2610
2611   TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2612   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2613   mult2->LoadHistograms("Multiplicity");
2614
2615   // "normal" result
2616   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2617
2618   if (chi2)
2619   {
2620     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2621     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2622   }
2623   else
2624     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2625
2626   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2627   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2628
2629   TH1* syst[2];
2630
2631   // change of pt spectrum (down)
2632   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2633   AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2634   mult3->LoadHistograms("Multiplicity");
2635   mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2636   if (chi2)
2637   {
2638     mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2639     mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2640   }
2641   else
2642     mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2643   syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2644
2645   // change of pt spectrum (up)
2646   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2647   AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2648   mult4->LoadHistograms("Multiplicity");
2649   mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2650   if (chi2)
2651   {
2652     mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2653     mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2654   }
2655   else
2656     mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2657   syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
2658
2659   DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
2660
2661   Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2662   Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2663 }
2664
2665 TH1* SystematicsSummary(Bool_t tpc = 0, Bool_t nsd = kTRUE)
2666 {
2667   Int_t nEffects = 7;
2668
2669   TH1* effects[10];
2670   const char** names = 0;
2671   Int_t colors[] = { 1, 2, 4, 1, 2, 4 };
2672   Int_t styles[] = { 1, 2, 3, 1, 2, 3 };
2673   Int_t widths[] = { 1, 1, 1, 2, 2, 2 };
2674   Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
2675
2676   TH1* dummy = new TH2F("dummy", Form(";%s;Uncertainty", GetMultLabel()), 202, -1.5, 200.5, 100, 0, 0.4);
2677   dummy->SetStats(0);
2678
2679   for (Int_t i=0; i<nEffects; ++i)
2680     effects[i] = new TH1F("SystematicsSummary", Form(";%s;Uncertainty", GetMultLabel()), 201, -0.5, 200.5);
2681
2682   if (tpc)
2683   {
2684     SetTPC();
2685
2686     const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2687     names = namesTPC;
2688
2689     // method
2690     TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2691     TH1* hist = (TH1*) file->Get("errorBoth");
2692
2693     // smooth a bit, but skip 0 bin
2694     effects[0]->SetBinContent(2, hist->GetBinContent(2));
2695     for (Int_t i=3; i<=200; ++i)
2696       effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2697
2698     // relative x-section
2699     effects[1]->SetBinContent(2, 0.005);
2700     effects[1]->SetBinContent(3, 0.0025);
2701     effects[1]->SetBinContent(4, 0.0025);
2702
2703     // particle composition
2704     for (Int_t i=2; i<=101; ++i)
2705     {
2706       if (i < 41)
2707       {
2708         effects[2]->SetBinContent(i, 0.01);
2709       }
2710       else if (i < 76)
2711       {
2712         effects[2]->SetBinContent(i, 0.02);
2713       }
2714       else
2715         effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2716     }
2717
2718     // pt cut off (only tpc)
2719     for (Int_t i=2; i<=101; ++i)
2720     {
2721       if (i < 11)
2722       {
2723         effects[3]->SetBinContent(i, 0.05);
2724       }
2725       else if (i < 51)
2726       {
2727         effects[3]->SetBinContent(i, 0.01);
2728       }
2729       else
2730         effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2731     }
2732
2733     // track selection (only tpc)
2734     for (Int_t i=2; i<=101; ++i)
2735       effects[4]->SetBinContent(i, 0.03);
2736
2737     // secondaries
2738     for (Int_t i=2; i<=101; ++i)
2739       effects[5]->SetBinContent(i, 0.01);
2740
2741     // pt spectrum
2742     for (Int_t i=2; i<=101; ++i)
2743     {
2744       if (i < 21)
2745       {
2746         effects[6]->SetBinContent(i, 0.05);
2747       }
2748       else if (i < 51)
2749       {
2750         effects[6]->SetBinContent(i, 0.02);
2751       }
2752       else
2753         effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2754     }
2755
2756   }
2757   else
2758   {
2759     nEffects = 5;
2760
2761     //const char* namesSPD[] = { "Particle composition",  "p_{t} cut-off", "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)"};
2762     const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)", "Particle composition",  "p_{t} cut-off"};
2763     names = namesSPD;
2764
2765     currentEffect = 0;
2766
2767     // method
2768     TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2769     TH1* hist = (TH1*) file->Get("errorBoth");
2770     hist->Scale(1.0 / sqrt(2));
2771
2772     // smooth a bit, but skip 0 bin
2773     /*effects[currentEffect]->SetBinContent(1, hist->GetBinContent(1));
2774     for (Int_t i=2; i<=201; ++i)
2775       effects[currentEffect]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);*/
2776     effects[currentEffect] = hist;
2777
2778     currentEffect++;
2779
2780     // relative x-section
2781     file = TFile::Open("uncertainty_xsection.root");
2782     effects[currentEffect++] = (TH1*) file->Get("effError_INEL");
2783     effects[currentEffect] = (TH1*) file->Get("effError_NSD");
2784     effects[currentEffect]->SetLineStyle(1);
2785     //effects[2]->SetBinContent(1, 0.20);
2786     //effects[2]->SetBinContent(2, 0.01);
2787     //effects[2]->SetBinContent(3, 0.002);
2788
2789     currentEffect++;
2790
2791     // particle composition
2792     effects[currentEffect]->SetBinContent(1, 0.16);
2793     for (Int_t i=2; i<=81; ++i)
2794     {
2795       effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * i / 81);
2796     }
2797
2798     currentEffect++;
2799
2800     // pt spectrum
2801     effects[currentEffect]->SetBinContent(1, 0.06);
2802     effects[currentEffect]->SetBinContent(2, 0.03);
2803     for (Int_t i=3; i<=81; ++i)
2804     {
2805       if (i <= 61)
2806       {
2807         effects[currentEffect]->SetBinContent(i, 0.01);
2808       }
2809       else if (i <= 81)
2810       {
2811         effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * (i - 61) / 20);
2812       }
2813     }
2814
2815 //     currentEffect++;
2816 //
2817 //     // material budget
2818 //     for (Int_t i=1; i<=81; ++i)
2819 //     {
2820 //       if (i < 5)
2821 //         effects[currentEffect]->SetBinContent(i, 0.05 - 0.01 * i);
2822 //       if (i > 51)
2823 //         effects[currentEffect]->SetBinContent(i, 0.05 * (i - 50) / 30);
2824 //     }
2825 //
2826     currentEffect++;
2827
2828   }
2829
2830   TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 500);
2831   canvas->SetRightMargin(0.05);
2832   canvas->SetTopMargin(0.05);
2833   //canvas->SetGridx();
2834   canvas->SetGridy();
2835   TLegend* legend = new TLegend(0.2, 0.4, 0.7, 0.4 + 0.5 * nEffects / 7);
2836   legend->SetFillColor(0);
2837   legend->SetTextSize(0.04);
2838   dummy->Draw();
2839   dummy->GetXaxis()->SetRangeUser(0, displayRange);
2840
2841   for (Int_t i=0; i<nEffects; ++i)
2842   {
2843     TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2844     /*current->Reset();
2845     for (Int_t j=0; j<nEffects-i; ++j)
2846       current->Add(effects[j]);*/
2847
2848     current->SetLineColor(colors[i]);
2849     current->SetLineStyle(styles[i]);
2850     current->SetLineWidth(widths[i]);
2851     //current->SetFillColor(colors[i]);
2852     current->SetMarkerColor(colors[i]);
2853     //current->SetMarkerStyle(markers[i]);
2854
2855     current->SetStats(kFALSE);
2856     current->GetYaxis()->SetRangeUser(0, 0.4);
2857     current->DrawCopy("SAME");
2858     legend->AddEntry(current, names[i]);
2859
2860     //TLatex* text = new TLatex(displayRange+2, current->GetBinContent(displayRange+1), names[i]);
2861     TLatex* text = new TLatex(displayRange+2, 0.1 - i * 0.02, names[i]);
2862     text->SetTextSize(0.04);
2863     text->SetTextColor(colors[i]);
2864     //text->Draw();
2865   }
2866
2867   // add total in square
2868   TH1* totalINEL = (TH1*) effects[0]->Clone("totalINEL");
2869   totalINEL->Reset();
2870   TH1* totalNSD = (TH1*) totalINEL->Clone("totalNSD");
2871
2872   for (Int_t i=0; i<nEffects; ++i)
2873   {
2874     //Printf("%d %f", i, effects[i]->GetBinContent(20));
2875     effects[i]->Multiply(effects[i]);
2876
2877     if (i != 2)
2878       totalINEL->Add(effects[i]);
2879     if (i != 1)
2880       totalNSD->Add(effects[i]);
2881   }
2882
2883   for (Int_t i=1; i<=totalINEL->GetNbinsX(); ++i)
2884   {
2885     totalINEL->SetBinError(i, 0);
2886     if (totalINEL->GetBinContent(i) > 0)
2887       totalINEL->SetBinContent(i, TMath::Min(sqrt(totalINEL->GetBinContent(i)), 1.0));
2888     totalNSD->SetBinError(i, 0);
2889     if (totalNSD->GetBinContent(i) > 0)
2890       totalNSD->SetBinContent(i, TMath::Min(sqrt(totalNSD->GetBinContent(i)), 1.0));
2891   }
2892
2893   //Printf("%f", total->GetBinContent(20));
2894
2895   totalINEL->SetMarkerStyle(5);
2896   totalINEL->SetMarkerColor(1);
2897   legend->AddEntry(totalINEL, "Total (INEL)", "P");
2898
2899   totalNSD->SetMarkerStyle(24);
2900   totalNSD->SetMarkerColor(2);
2901   legend->AddEntry(totalNSD, "Total (NSD)", "P");
2902
2903   Printf("total in bin 0 is INEL: %f NSD: %f", totalINEL->GetBinContent(1), totalNSD->GetBinContent(1));
2904   totalINEL->DrawCopy("SAME P"); //->SetBinContent(1, 0);
2905   totalNSD->DrawCopy("SAME P"); //->SetBinContent(1, 0);
2906
2907   legend->Draw();
2908
2909   canvas->SaveAs(canvas->GetName());
2910
2911   return (nsd) ? totalNSD : totalINEL;
2912 }
2913
2914 void finalPlot(Bool_t tpc = 0, Bool_t small = kFALSE)
2915 {
2916   loadlibs();
2917
2918   if (tpc)
2919     SetTPC();
2920
2921   //TH1* errorNSD = SystematicsSummary(tpc, 1);
2922
2923   TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", 800, 500);
2924   canvas->SetRightMargin(0.05);
2925   canvas->SetTopMargin(0.05);
2926   canvas->SetGridx();
2927   canvas->SetGridy();
2928
2929   legend = new TLegend(0.5, 0.6, 0.9, 0.8);
2930   legend->SetFillColor(0);
2931   legend->SetTextSize(0.04);
2932
2933   for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
2934   {
2935     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
2936     TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
2937     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2938
2939     DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
2940
2941     // normalize result
2942     //result->Scale(1.0 / result->Integral(2, displayRange));
2943
2944     result->GetXaxis()->SetRangeUser(0, displayRange);
2945     //result->SetBinContent(1, 0); result->SetBinError(1, 0);
2946     result->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (etaRange+1) * 0.5));
2947     result->SetMarkerStyle(0);
2948     result->SetLineColor(1);
2949     result->SetStats(kFALSE);
2950
2951     // systematic error
2952     TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
2953
2954     TH1* systError = (TH1*) result->Clone("systError");
2955     for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
2956       systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2957
2958     // change error drawing style
2959     systError->SetFillColor(15);
2960
2961     if (eventType == AliMultiplicityCorrection::kNSD)
2962     {
2963       result->SetLineColor(2);
2964       result->SetMarkerColor(2);
2965       result->SetMarkerStyle(5);
2966     }
2967
2968     canvas->cd();
2969     systError->DrawCopy(Form("E2 ][ %s", (eventType == AliMultiplicityCorrection::kINEL) ? "" : "SAME"));
2970     result->DrawCopy("SAME E ][");
2971     canvas->SetLogy();
2972
2973     legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic cross-section" : "NSD cross-section", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
2974   }
2975
2976   legend->Draw();
2977   /*
2978   //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
2979   TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
2980   text->SetFillColor(0);
2981   text->SetTextAlign(12);
2982   text->AddText("Systematic errors summed quadratically");
2983   text->AddText("0.6 million minimum bias events");
2984   text->AddText("corrected to inelastic events");
2985   text->Draw("B");
2986
2987   TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
2988   text2->SetFillColor(0);
2989   text2->SetTextAlign(12);
2990   text2->AddText("#sqrt{s} = 14 TeV");
2991   if (tpc)
2992   {
2993     text2->AddText("|#eta| < 0.9");
2994   }
2995   else
2996     text2->AddText("|#eta| < 2.0");
2997   text2->AddText("simulated data (PYTHIA)");
2998   text2->Draw("B");
2999
3000   if (tpc)
3001   {
3002     TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
3003     text3->SetNDC();
3004     text3->Draw();
3005   }
3006   else
3007   {
3008     TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
3009     text3->SetNDC();
3010     text3->Draw();
3011   }
3012
3013   // alice logo
3014   TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
3015   pad->Draw();
3016   pad->cd();
3017   TImage* img = TImage::Open("$HOME/alice.png");
3018   img->SetImageQuality(TAttImage::kImgBest);
3019   img->Draw();
3020
3021   canvas->Modified();
3022   */
3023
3024 /*  TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
3025   text->SetTextSize(0.04);
3026   text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
3027   text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
3028   text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
3029   text->Draw();*/
3030
3031
3032   canvas->SaveAs(canvas->GetName());
3033 }
3034
3035 void finalPlot2(Bool_t tpc = 0)
3036 {
3037   loadlibs();
3038
3039   if (tpc)
3040     SetTPC();
3041
3042   //TH1* errorNSD = SystematicsSummary(tpc, 1);
3043
3044   TCanvas* canvas = new TCanvas("finalPlot2.eps", "finalPlot2.eps", 800, 600);
3045   canvas->SetRightMargin(0.05);
3046   canvas->SetTopMargin(0.05);
3047   canvas->SetGridx();
3048   canvas->SetGridy();
3049   canvas->SetLogy();
3050
3051   legend = new TLegend(0.5, 0.7, 0.9, 0.9);
3052   legend->SetFillColor(0);
3053   legend->SetTextSize(0.04);
3054
3055   Int_t displayRanges[] = { 50, 80, 120 };
3056
3057   TH2* dummy = new TH2F("dummy", ";True multiplicity N_{ch};P(N_{ch})", 100, -0.5, displayRanges[2]+10, 1000, 5e-5, 5);
3058   dummy->SetStats(0);
3059   dummy->Draw();
3060
3061   TList list;
3062
3063   Int_t count = 0;
3064   for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
3065   {
3066     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
3067     //TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
3068     for (Int_t etaR = 0; etaR <= 2; etaR++)
3069     {
3070       TH1* result = mult->GetMultiplicityESDCorrected(etaR);
3071
3072       //DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
3073
3074       // normalize result
3075       result->Scale(TMath::Power(5, etaR) / result->Integral(1, displayRanges[etaR]));
3076
3077       result->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
3078       //result->SetBinContent(1, 0); result->SetBinError(1, 0);
3079       //result->SetTitle(Form(""));
3080       result->SetMarkerStyle(0);
3081       result->SetLineColor(1);
3082       result->SetLineWidth(2);
3083       //result->SetMarkerStyle(4);
3084       //result->SetStats(kFALSE);
3085
3086       // systematic error
3087       TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
3088
3089       TH1* systError = (TH1*) result->Clone("systError");
3090       // small hack until we have syst. errors for all eta ranges
3091       Float_t factor = 80.0 / displayRanges[etaR];
3092       for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
3093       {
3094         systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(factor * i));
3095       }
3096
3097       // change error drawing style
3098       systError->SetFillColor(15);
3099
3100       if (eventType == AliMultiplicityCorrection::kNSD)
3101       {
3102         result->SetLineColor(2);
3103         result->SetMarkerColor(2);
3104         result->SetMarkerStyle(5);
3105       }
3106
3107       canvas->cd();
3108       systError->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
3109       systError->DrawCopy("E2 ][ SAME");
3110       list.Add(result);
3111
3112       if (eventType == AliMultiplicityCorrection::kINEL)
3113       {
3114         TLatex* Tl = new TLatex;
3115         Tl->SetTextSize(0.04);
3116         //Tl->SetBit(TLatex::kTextNDC);
3117         Tl->SetTextAlign(32);
3118
3119         Float_t etaRangeArr[] = { 0.5, 1.0, 1.4 };
3120         TString tmpStr;
3121         tmpStr.Form("|#eta| < %.1f (x %d)", etaRangeArr[etaR], (Int_t) TMath::Power(5, etaR));
3122         if (etaR == 0)
3123           Tl->DrawLatex(36, result->GetBinContent(41), tmpStr);
3124         if (etaR == 1)
3125         {
3126           Tl->SetTextAlign(12);
3127           Tl->DrawLatex(82, result->GetBinContent(81), tmpStr);
3128         }
3129         if (etaR == 2)
3130         {
3131           Tl->SetTextAlign(12);
3132           Tl->DrawLatex(106, result->GetBinContent(101), tmpStr);
3133         }
3134       }
3135
3136       if (etaR == 0)
3137         legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic events" : "NSD events", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
3138
3139       count++;
3140     }
3141   }
3142
3143   for (Int_t i=0; i<list.GetEntries(); i++)
3144     ((TH1*) list.At(i))->DrawCopy("SAME E ][");
3145
3146   legend->Draw();
3147
3148   canvas->SaveAs(canvas->GetName());
3149 }
3150
3151 TMatrixD* NonInvertable()
3152 {
3153   const Int_t kSize = 5;
3154
3155   TMatrixD matrix(kSize, kSize);
3156   for (Int_t x=0; x<kSize; x++)
3157   {
3158     for (Int_t y=0; y<kSize; y++)
3159     {
3160       if (x == y)
3161       {
3162         if (x == 0 || x == kSize -1)
3163         {
3164           matrix(x, y) = 0.75;
3165         }
3166         else
3167           matrix(x, y) = 0.5;
3168       }
3169       else if (TMath::Abs(x - y) == 1)
3170       {
3171         matrix(x, y) = 0.25;
3172       }
3173     }
3174   }
3175
3176   matrix.Print();
3177
3178   //TMatrixD inverted(matrix);
3179   //inverted.Invert();
3180
3181   //inverted.Print();
3182
3183   return new TMatrixD(matrix);
3184 }
3185
3186 void BlobelUnfoldingExample()
3187 {
3188   const Int_t kSize = 20;
3189
3190   TMatrixD matrix(kSize, kSize);
3191   for (Int_t x=0; x<kSize; x++)
3192   {
3193     for (Int_t y=0; y<kSize; y++)
3194     {
3195       if (x == y)
3196       {
3197         if (x == 0 || x == kSize -1)
3198         {
3199           matrix(x, y) = 0.75;
3200         }
3201         else
3202           matrix(x, y) = 0.5;
3203       }
3204       else if (TMath::Abs(x - y) == 1)
3205       {
3206         matrix(x, y) = 0.25;
3207       }
3208     }
3209   }
3210
3211   //matrix.Print();
3212
3213   TMatrixD inverted(matrix);
3214   inverted.Invert();
3215
3216   //inverted.Print();
3217
3218   TH1F* inputDist = new TH1F("inputDist", ";t;Entries", kSize, -0.5, (Float_t) kSize - 0.5);
3219   TVectorD inputDistVector(kSize);
3220   TH1F* unfolded = inputDist->Clone("unfolded");
3221   TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
3222   measuredIdealDist->SetTitle(";m;Entries");
3223   TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
3224
3225   TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
3226   // norm: 1/(sqrt(2pi)sigma)
3227   gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
3228   //gaus->Print();
3229
3230   for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
3231   {
3232     Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
3233     inputDist->SetBinContent(x, value);
3234     inputDistVector(x-1) = value;
3235   }
3236
3237   TVectorD measuredDistIdealVector = matrix * inputDistVector;
3238
3239   for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
3240     measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
3241
3242   measuredDist->FillRandom(measuredIdealDist, 10000);
3243
3244   // fill error matrix before scaling
3245   TMatrixD covarianceMatrixMeasured(kSize, kSize);
3246   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3247     covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
3248
3249   TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
3250   //covarianceMatrix.Print();
3251
3252   TVectorD measuredDistVector(kSize);
3253   for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
3254     measuredDistVector(x-1) = measuredDist->GetBinContent(x);
3255
3256   TVectorD unfoldedVector = inverted * measuredDistVector;
3257   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3258     unfolded->SetBinContent(x, unfoldedVector(x-1));
3259
3260   TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1200, 600);
3261   canvas->SetTopMargin(0.05);
3262   canvas->Divide(2, 1);
3263
3264   canvas->cd(1);
3265   canvas->cd(1)->SetLeftMargin(0.15);
3266   canvas->cd(1)->SetRightMargin(0.05);
3267   canvas->cd(1)->SetTopMargin(0.05);
3268   gPad->SetGridx();
3269   gPad->SetGridy();
3270   measuredDist->GetYaxis()->SetRangeUser(-600, 2799);
3271   measuredDist->GetYaxis()->SetTitleOffset(1.9);
3272   measuredDist->SetStats(0);
3273   measuredDist->DrawCopy();
3274   gaus->Draw("SAME");
3275
3276   canvas->cd(2);
3277   canvas->cd(2)->SetLeftMargin(0.15);
3278   canvas->cd(2)->SetRightMargin(0.05);
3279   canvas->cd(2)->SetTopMargin(0.05);
3280   gPad->SetGridx();
3281   gPad->SetGridy();
3282   unfolded->GetYaxis()->SetRangeUser(-600, 2799);
3283   unfolded->GetYaxis()->SetTitleOffset(1.9);
3284   unfolded->SetStats(0);
3285   unfolded->DrawCopy();
3286   gaus->Draw("SAME");
3287
3288   canvas->SaveAs("BlobelUnfoldingExample.eps");
3289
3290   return;
3291
3292   // now unfold this with Bayesian
3293   loadlibs();
3294
3295   // fill a multiplicity object
3296   mult = new AliMultiplicityCorrection("mult", "mult");
3297   for (Int_t x=0; x<kSize; x++)
3298   {
3299     mult->GetMultiplicityVtx(0)->SetBinContent(1, x+1, inputDistVector(x));
3300     mult->GetMultiplicityESD(0)->SetBinContent(1, x+1, measuredDistVector(x)*10000);
3301     for (Int_t y=0; y<kSize; y++)
3302       mult->GetCorrelation(0)->SetBinContent(1, x+1, y+1, matrix(x, y));
3303   }
3304
3305   //mult->DrawHistograms();
3306
3307   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 0);
3308   //mult->SetCreateBigBin(kFALSE);
3309   mult->ApplyMinuitFit(0, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); //hist2->ProjectionY("mymchist"));
3310
3311   //mult->ApplyBayesianMethod(0, kFALSE, AliMultiplicityCorrection::kTrVtx, 0, -1, 0, kFALSE);
3312
3313   mult->DrawComparison("BlobelExample", 0, kFALSE, kTRUE, mult->GetMultiplicityVtx(0)->ProjectionY("mcmchist", 1, mult->GetMultiplicityVtx(0)->GetNbinsX()));
3314
3315
3316 }
3317
3318 void E735Fit()
3319 {
3320   TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
3321   fCurrentESD->Sumw2();
3322
3323   // Open the input stream
3324   ifstream in;
3325   in.open("e735data.txt");
3326
3327   while(in.good())
3328   {
3329     Float_t x, y, ye;
3330     in >> x >> y >> ye;
3331
3332     //Printf("%f %f %f", x, y, ye);
3333     fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
3334     fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
3335   }
3336
3337   in.close();
3338
3339   //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
3340
3341   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
3342
3343   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])");
3344   func->SetParNames("scaling", "averagen", "k");
3345   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
3346   func->SetParLimits(1, 0.001, 1000);
3347   func->SetParLimits(2, 0.001, 1000);
3348   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
3349
3350   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
3351   lognormal->SetParNames("scaling", "mean", "sigma");
3352   lognormal->SetParameters(1, 1, 1);
3353   lognormal->SetParLimits(0, 0, 10);
3354   lognormal->SetParLimits(1, 0, 100);
3355   lognormal->SetParLimits(2, 1e-3, 10);
3356
3357   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
3358   fCurrentESD->SetStats(kFALSE);
3359   fCurrentESD->SetMarkerStyle(0);
3360   fCurrentESD->SetLineColor(1);
3361   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
3362   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
3363   fCurrentESD->Draw("");
3364   fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
3365   fCurrentESD->Fit(func, "0", "", 0, 150);
3366   func->SetRange(0, 250);
3367   func->Draw("SAME");
3368   printf("chi2 = %f\n", func->GetChisquare());
3369
3370   fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
3371   lognormal->SetLineColor(2);
3372   lognormal->SetLineStyle(2);
3373   lognormal->SetRange(0, 250);
3374   lognormal->Draw("SAME");
3375
3376   gPad->SetLogy();
3377
3378   canvas->SaveAs("E735Fit.eps");
3379 }
3380
3381 void DifferentSamples()
3382 {
3383   loadlibs();
3384
3385   Int_t n = 2;
3386   const char* filesChi2[] = { "chi2_100k_1.root", "chi2_100k_2.root" };
3387   const char* filesBayesian[] = { "bayesian_100k_1.root", "bayesian_100k_2.root" };
3388
3389   TCanvas* canvas = new TCanvas("DifferentSamples", "DifferentSamples", 1200, 600);
3390   canvas->Divide(2, 1);
3391
3392   legend = new TLegend(0.15, 0.7, 0.65, 0.9);
3393   legend->SetFillColor(0);
3394   legend->SetTextSize(0.04);
3395
3396   for (Int_t i=0; i<n; i++)
3397   {
3398     AliMultiplicityCorrection* chi2 = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3399     TFile::Open(filesChi2[i]);
3400     chi2->LoadHistograms("Multiplicity");
3401
3402     AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3403     TFile::Open(filesBayesian[i]);
3404     bayesian->LoadHistograms("Multiplicity");
3405
3406     chi2Hist = chi2->GetMultiplicityESDCorrected(etaRange);
3407     bayesianHist = bayesian->GetMultiplicityESDCorrected(etaRange);
3408
3409     mc = chi2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
3410
3411     // normalize and divide
3412     chi2Hist->Scale(1.0 / chi2Hist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3413     bayesianHist->Scale(1.0 / bayesianHist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3414
3415     chi2Hist->Divide(mc, chi2Hist);
3416     bayesianHist->Divide(mc, bayesianHist);
3417
3418     canvas->cd(i+1);
3419     gPad->SetTopMargin(0.05);
3420     gPad->SetRightMargin(0.05);
3421     //gPad->SetLeftMargin(0.12);
3422     gPad->SetGridx();
3423     gPad->SetGridy();
3424
3425     chi2Hist->GetXaxis()->SetRangeUser(0, displayRange);
3426     chi2Hist->GetYaxis()->SetTitleOffset(1.3);
3427     chi2Hist->SetStats(0);
3428     chi2Hist->SetTitle(Form(";%s;MC / unfolded", GetMultLabel()));
3429     chi2Hist->GetYaxis()->SetRangeUser(0.2, 1.8);
3430     chi2Hist->Draw("HIST");
3431
3432     for (Int_t x=1; x<=bayesianHist->GetNbinsX(); x++)
3433       bayesianHist->SetBinError(x, 1e-6);
3434
3435     bayesianHist->SetLineColor(2);
3436     bayesianHist->SetMarkerColor(2);
3437     bayesianHist->SetMarkerStyle(5);
3438     bayesianHist->Draw("HIST E SAME");
3439
3440     if (i == 0)
3441     {
3442       legend->AddEntry(chi2Hist, "#chi^{2}-minimization", "L");
3443       legend->AddEntry(bayesianHist, "Bayesian unfolding", "LP");
3444     }
3445     legend->Draw();
3446   }
3447
3448   canvas->SaveAs("DifferentSamples.eps");
3449 }
3450
3451 void PileUp()
3452 {
3453   loadlibs();
3454
3455   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
3456   hist2d = mult->GetMultiplicityMC(etaRange, AliMultiplicityCorrection::kINEL);
3457   mult1 = hist2d->ProjectionY("mult1", 1, hist2d->GetNbinsX());
3458
3459   conv = (TH1*) mult1->Clone("conv");
3460   conv->Reset();
3461
3462   mult1->Scale(1.0 / mult1->Integral());
3463
3464   for (Int_t i=1; i<=mult1->GetNbinsX(); i++)
3465     for (Int_t j=1; j<=mult1->GetNbinsX(); j++)
3466       conv->Fill(mult1->GetBinCenter(i)+mult1->GetBinCenter(j), mult1->GetBinContent(i) * mult1->GetBinContent(j));
3467
3468   conv->Scale(1.0 / conv->Integral());
3469
3470   c = new TCanvas("c", "c", 800, 500);
3471   gPad->SetLogy();
3472   gPad->SetTopMargin(0.05);
3473   gPad->SetRightMargin(0.05);
3474   mult1->SetTitle(Form(";%s;Probability", GetMultLabel()));
3475   mult1->SetStats(0);
3476   gPad->SetGridx();
3477   gPad->SetGridy();
3478   mult1->Draw();
3479   mult1->GetYaxis()->SetRangeUser(1e-7, 2 * mult1->GetMaximum());
3480   mult1->GetXaxis()->SetRangeUser(0, displayRange);
3481   mult1->GetXaxis()->SetTitleOffset(1.15);
3482   conv->SetLineColor(2);
3483   conv->SetMarkerColor(2);
3484   conv->SetMarkerStyle(5);
3485   conv->DrawCopy("SAME P");
3486
3487   conv->Scale(0.00058);
3488   conv->DrawCopy("SAME P");
3489
3490   legend = new TLegend(0.73, 0.73, 0.93, 0.93);
3491   legend->SetFillColor(0);
3492   legend->SetTextSize(0.04);
3493   legend->AddEntry(mult1, "1 collision");
3494   legend->AddEntry(conv, "2 collisions", "P");
3495   legend->Draw();
3496
3497   c->SaveAs("pileup.eps");
3498
3499   new TCanvas;
3500   conv->Divide(mult1);
3501   conv->Draw();
3502 }
3503
3504 void TestErrorDetermination(Int_t nRandomizations)
3505 {
3506   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);
3507   func->SetParNames("scaling", "averagen", "k");
3508   func->SetParameters(1, 15, 2);
3509
3510   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);
3511   func2->SetParNames("scaling", "averagen", "k");
3512   func2->SetParLimits(0, 0.5, 2);
3513   func2->SetParLimits(1, 1, 50);
3514   func2->SetParLimits(2, 1, 10);
3515   func2->SetParameters(1, 15, 2);
3516   //func2->FixParameter(0, 1);
3517
3518   //new TCanvas; func->Draw("L");
3519
3520   hist1 = new TH1F("hist1", "", 100, 0.5, 100.5);
3521   hist2 = new TH1F("hist2", "", 100, 0.5, 100.5);
3522   hist1->Sumw2();
3523
3524   TH1* params[3];
3525   params[0] = new TH1F("param_0", Form("param_%d", 0), 100, 0.95, 1.05);
3526   params[1] = new TH1F("param_1", Form("param_%d", 1), 100, 14, 16);
3527   params[2] = new TH1F("param_2", Form("param_%d", 2), 100, 1.8, 2.2);
3528
3529   const Int_t nTries = 1000;
3530   for (Int_t i=0; i<nTries; i++)
3531   {
3532     hist1->Reset();
3533
3534     if (nRandomizations == 1)
3535     {
3536       hist1->FillRandom("nbd", 10000);
3537     }
3538     else if (nRandomizations == 2)
3539     {
3540       hist2->Reset();
3541       hist2->FillRandom("nbd", 10000);
3542       hist1->FillRandom(hist2, 10000);
3543     }
3544     else if (nRandomizations == 3)
3545     {
3546       hist2->Reset();
3547       hist1->FillRandom("nbd", 10000);
3548       hist2->FillRandom(hist1, 10000);
3549       hist1->Reset();
3550       hist1->FillRandom(hist2, 10000);
3551     }
3552     else
3553       return;
3554
3555     //new TCanvas; hist1->Draw();
3556
3557     hist1->Scale(1.0 / hist1->Integral());
3558     hist1->Fit(func2, "NQ");
3559     hist1->Fit(func2, "NQ");
3560     for (Int_t j=0; j<3; j++)
3561       params[j]->Fill(func2->GetParameter(j));
3562   }
3563
3564   for (Int_t j=0; j<3; j++)
3565   {
3566     new TCanvas; params[j]->Draw();
3567     params[j]->Fit("gaus");
3568     Printf("sigma of param %d if %f", j, ((TF1*) params[j]->FindObject("gaus"))->GetParameter(2));
3569   }
3570 }
3571
3572 void DrawRawDistributions(const char* fileName = "multiplicityESD.root")
3573 {
3574   loadlibs();
3575
3576   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileName);
3577
3578   c = new TCanvas("c", "c", 600, 600);
3579
3580   dummy = new TH2F("dummy", ";measured multiplicity", 100, -0.5, 149.5, 100, 0.5, 4e4);
3581   dummy->SetStats(0);
3582   dummy->Draw();
3583   gPad->SetGridx();
3584   gPad->SetGridy();
3585   gPad->SetLogy();
3586
3587   Int_t colors[] = { 1, 2, 4 };
3588
3589   for (Int_t i=2; i>=0; i--)
3590   {
3591     hist = mult->GetMultiplicityESD(i)->ProjectionY();
3592
3593     hist->SetLineColor(colors[i]);
3594     hist->DrawCopy("SAME");
3595   }
3596
3597
3598 }