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[] = { "LHC09b9_0.9TeV_0T/mb1/spd/multiplicity.root", "LHC09b8_0.9TeV_0.5T/mb1/tpc/multiplicity.root" };
1480   //const char* fileName[] = { "spd/multiplicity.root", "tpc/multiplicity.root" };
1481   //const char* fileName[] = { "multiplicity.root", "multiplicity.root" };
1482   Float_t etaRangeArr[] = {0.49, 0.9};
1483   const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1484
1485   TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1486   canvas->Divide(2, 1);
1487
1488   TCanvas* canvas3 = new TCanvas("EfficiencySpecies_comb", "EfficiencySpecies_comb", 600, 600);
1489   gPad->SetGridx();
1490   gPad->SetGridy();
1491   gPad->SetRightMargin(0.05);
1492   gPad->SetTopMargin(0.05);
1493
1494   TLegend* legends[2];
1495
1496   for (Int_t loop=0; loop<2; ++loop)
1497   {
1498     Printf("%s", fileName[loop]);
1499
1500     TCanvas* canvas2 = new TCanvas(Form("EfficiencySpecies_%d", loop), Form("EfficiencySpecies_%d", loop), 600, 600);
1501     gPad->SetGridx();
1502     gPad->SetGridy();
1503     gPad->SetRightMargin(0.05);
1504     gPad->SetTopMargin(0.05);
1505
1506     AliCorrection* correction[8];
1507
1508     canvas->cd(loop+1);
1509
1510     gPad->SetGridx();
1511     gPad->SetGridy();
1512     gPad->SetRightMargin(0.05);
1513
1514     TLegend* legend = new TLegend(0.6, 0.4, 0.85, 0.6);
1515     legend->SetFillColor(0);
1516     legend->SetEntrySeparation(0.2);
1517     legend->SetTextSize(gStyle->GetTextSize());
1518
1519     legends[loop] = new TLegend(0.4+loop*0.3, 0.2, 0.6+loop*0.3, 0.5);
1520     legends[loop]->SetFillColor(0);
1521     legends[loop]->SetEntrySeparation(0.2);
1522     legends[loop]->SetTextSize(gStyle->GetTextSize());
1523     legends[loop]->SetHeader((loop == 0) ? "SPD" : "TPC");
1524
1525     Float_t below = 0;
1526     Float_t total = 0;
1527
1528     TFile* file = TFile::Open(fileName[loop]);
1529     if (!file)
1530     {
1531       Printf("Could not open %s", fileName[loop]);
1532       return;
1533     }
1534
1535     Float_t sumGen = 0;
1536     Float_t sumMeas = 0;
1537
1538     Float_t sumGenAbove02 = 0;
1539     Float_t sumMeasAbove02 = 0;
1540
1541     for (Int_t i=0; i<3; ++i)
1542     {
1543       Printf("correction %d", i);
1544
1545       TString name; name.Form("correction_%d", i);
1546       correction[i] = new AliCorrection(name, name);
1547       correction[i]->LoadHistograms();
1548
1549       TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1550       TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1551
1552       // limit vtx axis
1553       Float_t vtxRange = 3.9;
1554       gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1555       meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1556
1557       // 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)
1558       /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1559         for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1560         {
1561           gene->SetBinContent(x, 0, z, 0);
1562           gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1563           meas->SetBinContent(x, 0, z, 0);
1564           meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1565         }*/
1566
1567       // limit eta axis
1568       Float_t etaBegin = -etaRangeArr[loop];
1569       Float_t etaEnd   = etaRangeArr[loop];
1570       //etaBegin = 0.01;
1571       //etaEnd = -0.01;
1572       gene->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
1573       meas->GetYaxis()->SetRangeUser(etaBegin, etaEnd);
1574
1575       TH1* genePt = gene->Project3D(Form("z_%d", i));
1576       TH1* measPt = meas->Project3D(Form("z_%d", i));
1577
1578       if (i == 2)
1579       {
1580         Printf("WARNING: Rebinning for protons!");
1581         genePt->Rebin(2);
1582         measPt->Rebin(2);
1583       }
1584
1585       genePt->Sumw2();
1586       measPt->Sumw2();
1587
1588       for (Int_t x=0; x<=genePt->GetNbinsX()+1; x++)
1589       {
1590         genePt->SetBinError(x, TMath::Sqrt(genePt->GetBinContent(x)));
1591         measPt->SetBinError(x, TMath::Sqrt(measPt->GetBinContent(x)));
1592       }
1593
1594       sumGen += genePt->Integral();
1595       sumMeas += measPt->Integral();
1596
1597       sumGenAbove02 += genePt->Integral(genePt->GetXaxis()->FindBin(0.21), genePt->GetNbinsX());
1598       sumMeasAbove02 += measPt->Integral(genePt->GetXaxis()->FindBin(0.21), genePt->GetNbinsX());
1599
1600       TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1601       effPt->Reset();
1602       effPt->Divide(measPt, genePt, 1, 1, "B");
1603
1604       Int_t bin = 0;
1605       for (bin=20; bin>=1; bin--)
1606       {
1607         if (effPt->GetBinContent(bin) < 0.5)
1608           break;
1609       }
1610
1611       Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1612
1613       Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1614       Printf("%.4f of the particles are below that momentum", fraction);
1615
1616       below += genePt->Integral(1, bin);
1617       total += genePt->Integral();
1618
1619       effPt->SetLineColor(color[i]);
1620       effPt->SetMarkerColor(color[i]);
1621       effPt->SetMarkerStyle(marker[i]);
1622       effPt->SetMarkerSize(2);
1623
1624       effPt->GetXaxis()->SetRangeUser(0, 1);
1625       effPt->GetYaxis()->SetRangeUser(0.001, 1);
1626
1627       effPt->GetXaxis()->SetTitleOffset(1.1);
1628       effPt->GetYaxis()->SetTitleOffset(1.2);
1629
1630       effPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1631
1632       effPt->SetStats(kFALSE);
1633       effPt->SetTitle(titles[loop]);
1634       effPt->GetYaxis()->SetTitle("Efficiency");
1635
1636       canvas->cd(loop+1);
1637       effPt->DrawCopy((i == 0) ? "" : "SAME");
1638
1639       canvas2->cd();
1640       effPt->SetTitle("");
1641       effPt->DrawCopy((i == 0) ? "" : "SAME");
1642
1643       canvas3->cd();
1644       effPtClone = (TH1*) effPt->Clone("effPtClone");
1645       effPtClone->SetMarkerStyle(marker[i]-4*loop);
1646       effPtClone->DrawCopy((i == 0 && loop == 0) ? "" : "SAME");
1647
1648       legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1649       legends[loop]->AddEntry(effPtClone, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")), "P");
1650       //legend2->AddEntry(effPt, Form("%s %s", (loop == 0) ? "SPD" : "TPC", ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}"))), "P");
1651
1652       if (addDecayStopped)
1653       {
1654         name.Form("correction_%d", i+4);
1655         corr = new AliCorrection(name, name);
1656         corr->LoadHistograms();
1657
1658         TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
1659         TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
1660
1661         // limit axes
1662         gene->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1663         meas->GetXaxis()->SetRangeUser(-vtxRange, vtxRange);
1664         gene->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1665         meas->GetYaxis()->SetRangeUser(-etaRangeArr[loop], etaRangeArr[loop]);
1666
1667         TH1* decayed = gene->Project3D(Form("z_%d", i+4));
1668         TH1* stopped = meas->Project3D(Form("z_%d", i+4));
1669
1670         Printf("%d: %d decayed, %d stopped, out of %d", i, (Int_t) decayed->Integral(), (Int_t) stopped->Integral(), (Int_t) genePt->Integral());
1671
1672         decayed->Divide(decayed, genePt, 1, 1, "B");
1673         stopped->Divide(stopped, genePt, 1, 1, "B");
1674
1675         decayed->SetMarkerStyle(20);
1676         stopped->SetMarkerStyle(21);
1677         stopped->SetMarkerColor(2);
1678
1679         new TCanvas(Form("all_%d_%d", loop, i), Form("all_%d_%d", loop, i), 600, 600);
1680         effPt->DrawCopy();
1681         decayed->DrawCopy("SAME");
1682         stopped->DrawCopy("SAME");
1683
1684         decayed->Add(stopped);
1685         decayed->Add(effPt);
1686         decayed->SetMarkerStyle(22);
1687         decayed->SetMarkerColor(4);
1688         decayed->DrawCopy("SAME");
1689       }
1690
1691     }
1692
1693     Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1694
1695     Printf("%f measured, %f generated, effiency: %f", sumGen, sumMeas, sumMeas / sumGen);
1696     Printf("Above 0.2 GeV/c: %f measured, %f generated, effiency: %f", sumGenAbove02, sumMeasAbove02, sumMeasAbove02 / sumGenAbove02);
1697
1698     canvas->cd(loop+1);
1699     legend->Draw();
1700
1701     canvas2->cd();
1702     legend->Draw();
1703     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
1704   }
1705
1706   canvas->cd();
1707   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1708
1709   canvas3->cd();
1710   legends[0]->Draw();
1711   legends[1]->Draw();
1712   canvas3->SaveAs(Form("%s.eps", canvas3->GetName()));
1713 }
1714
1715 void DrawpTCutOff()
1716 {
1717 /*
1718 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt_ref.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1719 mv unfolded.root chi2_ptref.root
1720 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1721 mv unfolded.root chi2_pt0.root
1722 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1723 mv unfolded.root chi2_pt1.root
1724 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt0_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1725 mv unfolded.root chi2_pt0_25.root
1726 aliroot -b -q $ALICE_ROOT/PWG0/multiplicity/correct.C'("multiplicityMC_pt1_25.root", "Multiplicity", "multiplicityESD.root", kTRUE)'
1727 mv unfolded.root chi2_pt1_25.root
1728 */
1729
1730   loadlibs();
1731
1732   TH1* results[10];
1733   const char* files[] = { "chi2_ptref.root", "chi2_pt0.root", "chi2_pt1.root", "chi2_pt0_25.root", "chi2_pt1_25.root"};
1734
1735   Int_t nMax = 5;
1736   for (Int_t i = 0; i<nMax; ++i)
1737   {
1738     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(files[i]);
1739     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1740   }
1741
1742   const char* legendStrings[] = { "Reduced 50%", "Enhanced 50%", "Reduced 25%", "Enhanced 25%" };
1743   DrawRatio(results[0], nMax-1, results+1, "LowMomentumSyst.eps", kFALSE, legendStrings);
1744 }
1745
1746 void ParticleSpeciesComparison()
1747 {
1748   loadlibs();
1749
1750   TH1* results[10];
1751   TH1* mc = 0;
1752
1753   // loop over cases (normal, enhanced/reduced ratios)
1754   Int_t nMax = 9;
1755   for (Int_t i = 0; i<nMax; ++i)
1756   {
1757     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(Form("chi2_species_%d.root", i), Form("Multiplicity_%d", i));
1758     if (i == 0)
1759       mc = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymchist", 1, 1);
1760     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1761   }
1762
1763   DrawResultRatio(mc, results[0], "ParticleSpeciesComparison1_1.eps");
1764
1765   for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1766   {
1767     results[0]->SetBinError(i, 0);
1768     mc->SetBinError(i, 0);
1769   }
1770
1771   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" };
1772
1773   DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1774
1775   //not valid: draw chi2 uncertainty on top!
1776   /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1777   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1778   errorHist->SetLineColor(1);
1779   errorHist->SetLineWidth(2);
1780   TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1781   for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1782   {
1783     errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1784     errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1785   }
1786
1787   errorHist->DrawCopy("SAME");
1788   errorHist2->DrawCopy("SAME");*/
1789
1790   //canvas->SaveAs(canvas->GetName());
1791
1792   DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1793
1794   //errorHist->DrawCopy("SAME");
1795   //errorHist2->DrawCopy("SAME");
1796
1797   //canvas2->SaveAs(canvas2->GetName());
1798 }
1799
1800 /*void ParticleSpeciesComparison2()
1801 {
1802   gSystem->Load("libPWG0base");
1803
1804   const char* fileNameMC = "multiplicityMC_400k_syst.root";
1805   const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1806   Bool_t chi2 = 0;
1807
1808   TFile::Open(fileNameMC);
1809   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1810   mult->LoadHistograms();
1811
1812   TH1* mc[10];
1813   TH1* results[10];
1814
1815   // loop over cases (normal, enhanced/reduced ratios)
1816   Int_t nMax = 7;
1817   for (Int_t i = 0; i<nMax; ++i)
1818   {
1819     TString folder;
1820     folder.Form("Multiplicity_%d", i);
1821
1822     AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1823
1824     TFile::Open(fileNameESD);
1825     mult2->LoadHistograms();
1826
1827     mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1828
1829     if (chi2)
1830     {
1831       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1832       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1833       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1834     }
1835     else
1836     {
1837       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1838       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1839     }
1840
1841     //Float_t averageRatio = 0;
1842     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1843
1844     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1845
1846     TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1847     mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1848
1849     //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1850     //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1851
1852     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1853   }
1854
1855   DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1856 }*/
1857
1858 TH1* Invert(TH1* eff)
1859 {
1860   // calculate corr = 1 / eff
1861
1862   TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1863   corr->Reset();
1864
1865   for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1866   {
1867     if (eff->GetBinContent(i) > 0)
1868     {
1869       corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1870       corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1871     }
1872   }
1873
1874   return corr;
1875 }
1876
1877 void TriggerVertexCorrection()
1878 {
1879   //
1880   // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1881   //
1882
1883   loadlibs();
1884
1885   TFile::Open(correctionFile);
1886   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1887   mult->LoadHistograms("Multiplicity");
1888
1889   TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1890   TH1* corrNSD = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kNSD));
1891   TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1892
1893   TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 500);
1894   gPad->SetGridx();
1895   gPad->SetGridy();
1896   gPad->SetTopMargin(0.05);
1897   gPad->SetRightMargin(0.05);
1898
1899   corrINEL->SetStats(kFALSE);
1900   corrINEL->GetXaxis()->SetRangeUser(0, 12);
1901   corrINEL->GetYaxis()->SetRangeUser(0, 8);
1902   corrINEL->SetTitle(Form(";%s;Correction factor", GetMultLabel()));
1903   corrINEL->SetMarkerStyle(22);
1904   corrINEL->Draw("PE");
1905
1906   corrMB->SetStats(kFALSE);
1907   corrMB->SetLineColor(2);
1908   corrMB->SetMarkerStyle(25);
1909   corrMB->SetMarkerColor(2);
1910   corrMB->Draw("SAME PE");
1911
1912   corrNSD->SetLineColor(4);
1913   corrNSD->SetMarkerStyle(24);
1914   corrNSD->SetMarkerColor(4);
1915   corrNSD->Draw("SAME PE");
1916
1917   Printf("       MB  INEL  NSD");
1918   Printf("bin 0: %f %f %f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1));
1919   Printf("bin 1: %f %f %f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2));
1920
1921   TLegend* legend = new TLegend(0.3, 0.6, 0.85, 0.85);
1922   legend->SetFillColor(0);
1923   legend->AddEntry(corrINEL, "Correction to inelastic sample");
1924   legend->AddEntry(corrNSD, "Correction to NSD sample");
1925   legend->AddEntry(corrMB, "Correction to triggered sample");
1926   legend->SetTextSize(0.04);
1927
1928   legend->Draw();
1929
1930   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1931 }
1932
1933 void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1934 {
1935   loadlibs();
1936
1937   TFile::Open(correctionFile);
1938   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1939   mult->LoadHistograms("Multiplicity");
1940
1941   TFile::Open(measuredFile);
1942   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1943   mult2->LoadHistograms("Multiplicity");
1944
1945   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1946
1947   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
1948
1949   AliUnfolding::SetNbins(70, 70);
1950   //AliUnfolding::SetNbins(35, 35);
1951   AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, 1e3);
1952   AliUnfolding::SetBayesianParameters(1, 10);
1953
1954   TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1955
1956   new TCanvas; errorMeasured->Draw();
1957   new TCanvas;
1958
1959   mult->GetMultiplicityESDCorrected(etaRange)->Scale(1.0 / mult->GetMultiplicityESDCorrected(etaRange)->Integral());
1960   mult->GetMultiplicityESDCorrected(etaRange)->Draw();
1961   mcHist->Scale(1.0 / mcHist->Integral());
1962   mcHist->SetLineColor(2);
1963   mcHist->Draw("SAME");
1964   gPad->SetLogy();
1965
1966   return;
1967
1968   TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1969
1970   TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1971
1972   if (!mc)
1973   {
1974     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1975     DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1976   }
1977
1978   TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
1979   errorResponse->Write();
1980   errorMeasured->Write();
1981   errorBoth->Write();
1982   file->Close();
1983 }
1984
1985 void DrawStatisticalUncertainty()
1986 {
1987   TFile::Open("StatisticalUncertainty.root");
1988
1989   errorResponse = (TH1*) gFile->Get("errorResponse");
1990   errorMeasured = (TH1*) gFile->Get("errorMeasured");
1991   errorBoth = (TH1*) gFile->Get("errorBoth");
1992
1993   TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1994   canvas->SetGridx();
1995   canvas->SetGridy();
1996   canvas->SetRightMargin(0.05);
1997   canvas->SetTopMargin(0.05);
1998
1999   errorResponse->SetLineColor(1);
2000   errorResponse->GetXaxis()->SetRangeUser(0, longDisplayRange);
2001   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
2002   errorResponse->SetStats(kFALSE);
2003   errorResponse->SetTitle(";true multiplicity;Uncertainty");
2004
2005   errorResponse->Draw();
2006
2007   errorMeasured->SetLineColor(2);
2008   errorMeasured->Draw("SAME");
2009
2010   errorBoth->SetLineColor(4);
2011   errorBoth->Draw("SAME");
2012
2013   Printf("Average errorResponse: %f", errorResponse->Integral(2, displayRange) / (displayRange - 1));
2014   Printf("Average errorMeasured: %f", errorMeasured->Integral(2, displayRange) /  (displayRange - 1));
2015   Printf("Average errorBoth: %f", errorBoth->Integral(2, displayRange) /  (displayRange - 1));
2016
2017   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2018 }
2019
2020 void StatisticalUncertaintyCompare(const char* det = "SPD")
2021 {
2022   TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
2023   TH1* errorResponse = (TH1*) file1->Get("errorResponse");
2024   TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
2025   TH1* errorBoth = (TH1*) file1->Get("errorBoth");
2026
2027   TString str;
2028   str.Form("StatisticalUncertaintyCompare%s", det);
2029
2030   TCanvas* canvas = new TCanvas(str, str, 800, 500);
2031   canvas->SetGridx();
2032   canvas->SetGridy();
2033   canvas->SetRightMargin(0.05);
2034   canvas->SetTopMargin(0.05);
2035
2036   errorResponse->Scale(1.0 / sqrt(2));
2037   errorMeasured->Scale(1.0 / sqrt(2));
2038   errorBoth->Scale(1.0 / sqrt(2));
2039
2040   errorResponse->SetLineColor(1);
2041   errorResponse->GetXaxis()->SetRangeUser(0, displayRange);
2042   errorResponse->GetYaxis()->SetRangeUser(0, 0.18);
2043   errorResponse->SetStats(kFALSE);
2044   errorResponse->GetYaxis()->SetTitleOffset(1.2);
2045   errorResponse->SetTitle(Form(";%s;#sqrt{2}^{-1} #sigma(unfolded - unfolded_{0}) / unfolded_{0}", GetMultLabel()));
2046
2047   errorResponse->Draw();
2048
2049   errorMeasured->SetLineColor(2);
2050   errorMeasured->Draw("SAME");
2051
2052   errorBoth->SetLineColor(4);
2053   errorBoth->Draw("SAME");
2054
2055   TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
2056   TH1* errorResponse2 = (TH1*) file2->Get("errorResponse");
2057   TH1* errorMeasured2 = (TH1*) file2->Get("errorMeasured");
2058   TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
2059
2060   errorResponse2->Scale(1.0 / sqrt(2));
2061   errorMeasured2->Scale(1.0 / sqrt(2));
2062   errorBoth2->Scale(1.0 / sqrt(2));
2063
2064   errorResponse2->SetLineStyle(2);
2065   errorResponse2->Draw("SAME");
2066
2067   errorMeasured2->SetLineColor(2);
2068   errorMeasured2->SetLineStyle(2);
2069   errorMeasured2->Draw("SAME");
2070
2071   errorBoth2->SetLineColor(4);
2072   errorBoth2->SetLineStyle(2);
2073   errorBoth2->Draw("SAME");
2074
2075   TLegend* legend = new TLegend(0.2, 0.5, 0.8, 0.9);
2076   legend->SetFillColor(0);
2077   legend->SetTextSize(0.04);
2078   legend->AddEntry(errorBoth, "Both (Bayesian unfolding)");
2079   legend->AddEntry(errorMeasured, "Measured (Bayesian unfolding)");
2080   legend->AddEntry(errorResponse, "Response matrix (Bayesian unfolding)");
2081   legend->AddEntry(errorBoth2, "Both (#chi^{2}-minimization)");
2082   legend->AddEntry(errorMeasured2, "Measured (#chi^{2}-minimization)");
2083   legend->AddEntry(errorResponse2, "Response matrix (#chi^{2}-minimization)");
2084   legend->Draw();
2085
2086   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2087 }
2088
2089 void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
2090 {
2091  const char* files[] = { "multiplicityMC_nd.root", "multiplicityMC_sd.root", "multiplicityMC_dd.root", "multiplicityMC_xsection.root" };
2092
2093   loadlibs();
2094
2095   TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 600);
2096   canvas->SetGridx();
2097   canvas->SetGridy();
2098   canvas->SetRightMargin(0.05);
2099   canvas->SetTopMargin(0.05);
2100
2101   AliMultiplicityCorrection* data[4];
2102   TH1* effArray[4];
2103   TH1* effErrorArray[2];
2104
2105   Int_t markers[] = { 24, 25, 26, 5 };
2106   //Int_t markers[] = { 2, 25, 24, 5 };
2107   Int_t colors[] = { 1, 2, 4, 6 };
2108   //Int_t colors[] = { 1, 1, 1, 1 };
2109
2110   //TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
2111   TLegend* legend = new TLegend(0.3, 0.3, 0.9, 0.6);
2112   legend->SetTextSize(0.04);
2113   legend->SetFillColor(0);
2114
2115   for (Int_t i=0; i<4; ++i)
2116   {
2117     TString name;
2118     name.Form("Multiplicity_%d", i);
2119
2120     TFile::Open(files[i]);
2121     data[i] = new AliMultiplicityCorrection(name, name);
2122
2123     if (i < 3)
2124     {
2125       data[i]->LoadHistograms("Multiplicity");
2126     }
2127     else
2128       data[i]->LoadHistograms("Multiplicity_0");
2129
2130     TH1* eff = 0;
2131     if (eventType == -1)
2132     {
2133       eff = (TH1*) data[i]->GetTriggerEfficiency(etaRange)->Clone(Form("eff_%d", i));
2134     }
2135     else
2136       eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
2137     effArray[i] = eff;
2138
2139     eff->GetXaxis()->SetRangeUser(0, 15);
2140     eff->GetYaxis()->SetRangeUser(0, 1.19);
2141     eff->SetStats(kFALSE);
2142     eff->GetXaxis()->SetTitle(GetMultLabel());
2143     eff->GetYaxis()->SetTitle("Efficiency");
2144     eff->SetTitle("");
2145     eff->SetLineColor(colors[i]);
2146     eff->SetMarkerColor(colors[i]);
2147     eff->SetMarkerStyle(markers[i]);
2148
2149     if (i == 3)
2150     {
2151       // once for INEL, once for NSD
2152       for (AliMultiplicityCorrection::EventType eventType2 = AliMultiplicityCorrection::kINEL; eventType2 <= AliMultiplicityCorrection::kNSD; eventType2++)
2153       {
2154         effDiff = (TH1*) data[i]->GetEfficiency(etaRange, eventType2)->Clone(Form("effDiff_%d", i));
2155
2156         for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2157           effDiff->SetBinError(bin, 0);
2158
2159         // loop over cross section combinations
2160         for (Int_t j=1; j<7; ++j)
2161         {
2162           AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
2163           mult->LoadHistograms(Form("Multiplicity_%d", j));
2164
2165           TH1* eff2 = mult->GetEfficiency(etaRange, eventType2);
2166
2167           for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2168           {
2169             // TODO we could also do asymmetric errors here
2170             Float_t deviation = TMath::Abs(effDiff->GetBinContent(bin) - eff2->GetBinContent(bin));
2171
2172             effDiff->SetBinError(bin, TMath::Max(effDiff->GetBinError(bin), (Double_t) deviation));
2173           }
2174         }
2175
2176         for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2177         {
2178           //if (eventType2 == AliMultiplicityCorrection::kINEL)
2179             //eff->SetBinError(bin, 0);
2180             //eff->SetBinError(bin, effDiff->GetBinError(bin));
2181           if (bin < 20 && effDiff->GetBinContent(bin) > 0)
2182             Printf("Bin %d: Error: %.2f", bin, 100.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2183         }
2184
2185         if (uncertainty) {
2186                 TH1* effError = (TH1*) effDiff->Clone(Form("effError_%s", (eventType2 == AliMultiplicityCorrection::kINEL) ? "INEL" : "NSD"));
2187           effErrorArray[eventType2 - AliMultiplicityCorrection::kINEL] = effError;
2188                 effError->Reset();
2189
2190                 for (Int_t bin=1; bin<=effDiff->GetNbinsX(); bin++)
2191                   if (effDiff->GetBinContent(bin) > 0)
2192                     effError->SetBinContent(bin, 1.0 * effDiff->GetBinError(bin) / effDiff->GetBinContent(bin));
2193
2194                 effError->SetLineColor(1);
2195                 effError->SetMarkerStyle(1);
2196
2197           if (eventType2 == AliMultiplicityCorrection::kNSD)
2198             effError->SetLineStyle(2);
2199
2200           effError->DrawCopy("SAME HIST");
2201         }
2202       }
2203     }
2204
2205     canvas->cd();
2206     if (i == 0)
2207     {
2208       eff->DrawCopy("P");
2209     }
2210     else
2211       eff->DrawCopy("SAME P");
2212
2213     legend->AddEntry(eff, (((i == 0) ? "Non-diffractive" : ((i == 1) ? "Single-diffractive" : ((i == 2) ? "Double-diffractive" : "Pythia combined")))));
2214   }
2215
2216   if (uncertainty)
2217   {
2218     legend->AddEntry(effErrorArray[0], "Relative syst. uncertainty: inelastic");
2219     legend->AddEntry(effErrorArray[1], "Relative syst. uncertainty: NSD");
2220
2221     file = TFile::Open("uncertainty_xsection.root", "RECREATE");
2222     effErrorArray[0]->Write();
2223     effErrorArray[1]->Write();
2224     file->Close();
2225   }
2226
2227   legend->Draw();
2228
2229   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2230 }
2231
2232 void ModelDependencyPlot()
2233 {
2234   loadlibs();
2235
2236   TFile::Open("multiplicityMC.root");
2237   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2238   mult->LoadHistograms("Multiplicity");
2239
2240   hist = mult->GetCorrelation(etaRange);
2241
2242   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
2243   {
2244     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
2245     {
2246       hist->SetBinContent(0, y, z, 0);
2247       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2248     }
2249   }
2250
2251   TH2* proj = (TH2*) hist->Project3D("zy");
2252
2253   TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 1200, 600);
2254
2255   canvas->Divide(2, 1);
2256
2257   canvas->cd(2);
2258   gPad->SetLogy();
2259   gPad->SetGridx();
2260   gPad->SetGridy();
2261   gPad->SetTopMargin(0.05);
2262   gPad->SetRightMargin(0.05);
2263
2264   Int_t selectedMult = 30;
2265   Int_t yMax = 9e4;
2266
2267   TH1* full = proj->ProjectionX("full");
2268   TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
2269
2270   full->SetStats(kFALSE);
2271   full->GetXaxis()->SetRangeUser(0, displayRange);
2272   full->GetYaxis()->SetRangeUser(5, yMax);
2273   full->SetTitle(";Multiplicity;Entries");
2274
2275   selected->SetLineColor(0);
2276   selected->SetMarkerColor(2);
2277   selected->SetMarkerStyle(5);
2278
2279   full->Draw();
2280   selected->Draw("SAME P");
2281
2282   TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
2283   legend->SetFillColor(0);
2284   legend->AddEntry(full, "True");
2285   legend->AddEntry(selected, "Measured");
2286   legend->Draw();
2287
2288   TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
2289   line->SetLineWidth(2);
2290   line->Draw();
2291
2292   canvas->cd(1);
2293   gPad->SetLogy();
2294   gPad->SetGridx();
2295   gPad->SetGridy();
2296   gPad->SetTopMargin(0.05);
2297   gPad->SetRightMargin(0.05);
2298
2299   full = proj->ProjectionY("full2");
2300   selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
2301
2302   full->SetStats(kFALSE);
2303   full->GetXaxis()->SetRangeUser(0, displayRange);
2304   full->GetYaxis()->SetRangeUser(5, yMax);
2305   full->SetTitle(";Multiplicity;Entries");
2306
2307   full->SetLineColor(0);
2308   full->SetMarkerColor(2);
2309   full->SetMarkerStyle(5);
2310
2311   full->Draw("P");
2312   selected->Draw("SAME");
2313
2314   legend->Draw();
2315
2316   line = new TLine(selectedMult, 5, selectedMult, yMax);
2317   line->SetLineWidth(2);
2318   line->Draw();
2319
2320   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2321 }
2322
2323 void SystematicpTSpectrum()
2324 {
2325   gSystem->Load("libPWG0base");
2326
2327   TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
2328   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2329   mult->LoadHistograms("Multiplicity");
2330
2331   TFile::Open("multiplicityMC_100k_syst.root");
2332   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2333   mult2->LoadHistograms("Multiplicity");
2334
2335   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2336   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2337   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2338
2339   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2340   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2341
2342   DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
2343 }
2344
2345 void FitPt(const char* fileName)
2346 {
2347   // needs a MC file from the dNdEta analysis
2348
2349   TFile::Open(fileName);
2350
2351   TH1* genePt = (TH1*) gFile->Get("fdNdpT");
2352
2353   genePt->SetTitle(";p_{T} (GeV/c);1/p_{T} dN_{ch}/dp_{T} (GeV/c)^{-2}");
2354   // number of events
2355   genePt->Scale(1.0 / 287800);
2356   // bin width
2357   genePt->Scale(1.0 / genePt->GetXaxis()->GetBinWidth(1));
2358
2359   genePt->GetXaxis()->SetRangeUser(0, 0.4);
2360
2361   TF1* func = new TF1("func", "[1]*x*exp(x*[0])");
2362   func->SetParameters(-1, 1);
2363
2364   genePt->SetMarkerStyle(25);
2365   genePt->SetTitle("");
2366   genePt->SetStats(kFALSE);
2367
2368   new TCanvas;
2369   genePt->DrawCopy("P");
2370   func->DrawCopy("SAME");
2371   gPad->SetLogy();
2372
2373   genePt->Fit(func, "0", "", 0, 0.25);
2374   genePt->Fit(func, "0", "", 0, 0.25);
2375
2376   TCanvas* canvas = new TCanvas("FitPt", "FitPt", 600, 600);
2377
2378   gPad->SetGridx();
2379   gPad->SetGridy();
2380   gPad->SetLeftMargin(0.13);
2381   gPad->SetRightMargin(0.05);
2382   gPad->SetTopMargin(0.05);
2383
2384   //genePt->GetXaxis()->SetRangeUser(0, 1);
2385   genePt->GetYaxis()->SetRangeUser(2, 200);
2386   genePt->GetYaxis()->SetTitleOffset(1.4);
2387   genePt->GetXaxis()->SetTitleOffset(1.1);
2388   genePt->DrawCopy("P");
2389   //func->SetRange(0, 0.3);
2390   func->DrawCopy("SAME");
2391   gPad->SetLogy();
2392
2393   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2394
2395   TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2396
2397   TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2398
2399   TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2400
2401   for (Int_t param=0; param<2; param++)
2402   {
2403     for (Int_t sign=0; sign<2; sign++)
2404     {
2405       TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign));  // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2406       func2->SetParameters(func->GetParameters());
2407       //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2408
2409       Float_t factor = ((sign == 0) ? 0.75 : 1.25);
2410       func2->SetParameter(param, func2->GetParameter(param) * factor);
2411       //func2->Print();
2412
2413       canvas->cd();
2414       func2->SetLineWidth(2);
2415       func2->SetLineColor(2);
2416       func2->SetLineStyle(2);
2417       func2->DrawCopy("SAME");
2418
2419       canvas2->cd();
2420       TH1* second = func2->GetHistogram();
2421       second->Divide(first);
2422       second->SetLineColor(param + 1);
2423       // set to 1 above 0.2 GeV
2424       for (Int_t bin=second->GetXaxis()->FindBin(0.20001); bin<=second->GetNbinsX(); bin++)
2425         second->SetBinContent(bin, 1);
2426       second->GetYaxis()->SetRangeUser(0, 2);
2427       second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2428       second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2429     }
2430   }
2431
2432   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2433   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2434
2435   file->Close();
2436 }
2437
2438 void DrawSystematicpT()
2439 {
2440   TFile* file = TFile::Open("SystematicpT.root");
2441
2442   TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2443   TH1* result2 = (TH1*) file->Get("result_unity");
2444
2445   TH1* mcHist[12];
2446   TH1* result[12];
2447
2448   Int_t nParams = 5;
2449
2450   for (Int_t id=0; id<nParams*2; ++id)
2451   {
2452     mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2453     result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2454   }
2455
2456   DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2457
2458   //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2459
2460   DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2461
2462   //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2463
2464   // does not make sense: mc is different
2465   //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2466 }
2467
2468 void SystematicpT(Bool_t chi2 = 1)
2469 {
2470   gSystem->Load("libPWG0base");
2471
2472   TFile::Open("ptspectrum900.root");
2473   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2474   mult->LoadHistograms("Multiplicity");
2475
2476   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2477
2478   TH1* mcHist[12];
2479   TH1* result[12];
2480
2481   Int_t nParams = 5;
2482
2483   for (Int_t param=0; param<nParams; param++)
2484   {
2485     for (Int_t sign=0; sign<2; sign++)
2486     {
2487       // calculate result with systematic effect
2488       TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2489       mult2->LoadHistograms("Multiplicity");
2490
2491       mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2492
2493       if (chi2)
2494       {
2495         mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2496         mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2497       }
2498       else
2499         mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2500
2501       Int_t id = param * 2 + sign;
2502
2503       mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2504       result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2505
2506       TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2507       DrawResultRatio(mcHist[id], result[id], tmp);
2508     }
2509   }
2510
2511   // calculate normal result
2512   TFile::Open("ptspectrum100_1.root");
2513   mult2->LoadHistograms("Multiplicity");
2514   TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2515
2516   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2517
2518   if (chi2)
2519   {
2520     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2521   }
2522   else
2523     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2524
2525   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2526
2527   TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2528   mcHist2->Write();
2529   result2->Write();
2530   for (Int_t id=0; id<nParams*2; ++id)
2531   {
2532     mcHist[id]->Write();
2533     result[id]->Write();
2534   }
2535   file->Close();
2536
2537   DrawSystematicpT();
2538 }
2539
2540 void DrawSystematicpT2()
2541 {
2542   //displayRange = 200;
2543
2544   // read from file
2545   TFile* file = TFile::Open("SystematicpT2.root");
2546   TH1* mcHist = (TH1*) file->Get("mymc");
2547   TH1* result[12];
2548   result[0] = (TH1*) file->Get("result_unity");
2549   Int_t nParams = 5;
2550   for (Int_t id=0; id<nParams*2; ++id)
2551     result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2552
2553   DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2554   DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2555   DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2556 }
2557
2558 void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2559 {
2560   gSystem->Load("libPWG0base");
2561
2562   if (tpc)
2563   {
2564     SetTPC();
2565     TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2566   }
2567   else
2568     TFile::Open("ptspectrum100_1.root");
2569
2570   AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2571   measured->LoadHistograms("Multiplicity");
2572   TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2573
2574   TH1* result[12];
2575
2576   Int_t nParams = 5;
2577
2578   // -1 = unity change, 0...4 parameters
2579   for (Int_t id=-1; id<nParams*2; id++)
2580   {
2581     Int_t param = id / 2;
2582     Int_t sign = id % 2;
2583
2584     TString idStr;
2585     if (id == -1)
2586     {
2587       idStr = "unity";
2588     }
2589     else
2590       idStr.Form("%d_%d", param, sign);
2591
2592     // calculate result with systematic effect
2593     if (tpc)
2594     {
2595       TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2596     }
2597     else
2598       TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2599
2600     AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2601     response->LoadHistograms("Multiplicity");
2602
2603     response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2604
2605     if (chi2)
2606     {
2607       response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2608       response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2609     }
2610     else
2611       response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2612
2613     result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2614
2615     TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2616     DrawResultRatio(mcHist, result[id+1], tmp);
2617   }
2618
2619   TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2620   mcHist->Write();
2621   for (Int_t id=0; id<nParams*2+1; ++id)
2622     result[id]->Write();
2623   file->Close();
2624
2625   DrawSystematicpT2();
2626 }
2627
2628 void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2629 {
2630   // only needed for TPC
2631   SetTPC();
2632
2633   gSystem->Load("libPWG0base");
2634
2635   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2636   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2637   mult->LoadHistograms("Multiplicity");
2638
2639   TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2640   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2641   mult2->LoadHistograms("Multiplicity");
2642
2643   // "normal" result
2644   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2645
2646   if (chi2)
2647   {
2648     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2649     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2650   }
2651   else
2652     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2653
2654   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2655   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2656
2657   TH1* syst[2];
2658
2659   // change of pt spectrum (down)
2660   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2661   AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2662   mult3->LoadHistograms("Multiplicity");
2663   mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2664   if (chi2)
2665   {
2666     mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2667     mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2668   }
2669   else
2670     mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2671   syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2672
2673   // change of pt spectrum (up)
2674   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2675   AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2676   mult4->LoadHistograms("Multiplicity");
2677   mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2678   if (chi2)
2679   {
2680     mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2681     mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2682   }
2683   else
2684     mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2685   syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
2686
2687   DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
2688
2689   Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2690   Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2691 }
2692
2693 TH1* SystematicsSummary(Bool_t tpc = 0, Bool_t nsd = kTRUE)
2694 {
2695   Int_t nEffects = 7;
2696
2697   TH1* effects[10];
2698   const char** names = 0;
2699   Int_t colors[] = { 1, 2, 4, 1, 2, 4 };
2700   Int_t styles[] = { 1, 2, 3, 1, 2, 3 };
2701   Int_t widths[] = { 1, 1, 1, 2, 2, 2 };
2702   Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
2703
2704   TH1* dummy = new TH2F("dummy", Form(";%s;Uncertainty", GetMultLabel()), 202, -1.5, 200.5, 100, 0, 0.4);
2705   dummy->SetStats(0);
2706
2707   for (Int_t i=0; i<nEffects; ++i)
2708     effects[i] = new TH1F("SystematicsSummary", Form(";%s;Uncertainty", GetMultLabel()), 201, -0.5, 200.5);
2709
2710   if (tpc)
2711   {
2712     SetTPC();
2713
2714     const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2715     names = namesTPC;
2716
2717     // method
2718     TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2719     TH1* hist = (TH1*) file->Get("errorBoth");
2720
2721     // smooth a bit, but skip 0 bin
2722     effects[0]->SetBinContent(2, hist->GetBinContent(2));
2723     for (Int_t i=3; i<=200; ++i)
2724       effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2725
2726     // relative x-section
2727     effects[1]->SetBinContent(2, 0.005);
2728     effects[1]->SetBinContent(3, 0.0025);
2729     effects[1]->SetBinContent(4, 0.0025);
2730
2731     // particle composition
2732     for (Int_t i=2; i<=101; ++i)
2733     {
2734       if (i < 41)
2735       {
2736         effects[2]->SetBinContent(i, 0.01);
2737       }
2738       else if (i < 76)
2739       {
2740         effects[2]->SetBinContent(i, 0.02);
2741       }
2742       else
2743         effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2744     }
2745
2746     // pt cut off (only tpc)
2747     for (Int_t i=2; i<=101; ++i)
2748     {
2749       if (i < 11)
2750       {
2751         effects[3]->SetBinContent(i, 0.05);
2752       }
2753       else if (i < 51)
2754       {
2755         effects[3]->SetBinContent(i, 0.01);
2756       }
2757       else
2758         effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2759     }
2760
2761     // track selection (only tpc)
2762     for (Int_t i=2; i<=101; ++i)
2763       effects[4]->SetBinContent(i, 0.03);
2764
2765     // secondaries
2766     for (Int_t i=2; i<=101; ++i)
2767       effects[5]->SetBinContent(i, 0.01);
2768
2769     // pt spectrum
2770     for (Int_t i=2; i<=101; ++i)
2771     {
2772       if (i < 21)
2773       {
2774         effects[6]->SetBinContent(i, 0.05);
2775       }
2776       else if (i < 51)
2777       {
2778         effects[6]->SetBinContent(i, 0.02);
2779       }
2780       else
2781         effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2782     }
2783
2784   }
2785   else
2786   {
2787     nEffects = 5;
2788
2789     //const char* namesSPD[] = { "Particle composition",  "p_{t} cut-off", "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)"};
2790     const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Relative cross-sections (INEL)", "Relative cross-sections (NSD)", "Particle composition",  "p_{t} cut-off"};
2791     names = namesSPD;
2792
2793     currentEffect = 0;
2794
2795     // method
2796     TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2797     TH1* hist = (TH1*) file->Get("errorBoth");
2798     hist->Scale(1.0 / sqrt(2));
2799
2800     // smooth a bit, but skip 0 bin
2801     /*effects[currentEffect]->SetBinContent(1, hist->GetBinContent(1));
2802     for (Int_t i=2; i<=201; ++i)
2803       effects[currentEffect]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);*/
2804     effects[currentEffect] = hist;
2805
2806     currentEffect++;
2807
2808     // relative x-section
2809     file = TFile::Open("uncertainty_xsection.root");
2810     effects[currentEffect++] = (TH1*) file->Get("effError_INEL");
2811     effects[currentEffect] = (TH1*) file->Get("effError_NSD");
2812     effects[currentEffect]->SetLineStyle(1);
2813     //effects[2]->SetBinContent(1, 0.20);
2814     //effects[2]->SetBinContent(2, 0.01);
2815     //effects[2]->SetBinContent(3, 0.002);
2816
2817     currentEffect++;
2818
2819     // particle composition
2820     effects[currentEffect]->SetBinContent(1, 0.16);
2821     for (Int_t i=2; i<=81; ++i)
2822     {
2823       effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * i / 81);
2824     }
2825
2826     currentEffect++;
2827
2828     // pt spectrum
2829     effects[currentEffect]->SetBinContent(1, 0.06);
2830     effects[currentEffect]->SetBinContent(2, 0.03);
2831     for (Int_t i=3; i<=81; ++i)
2832     {
2833       if (i <= 61)
2834       {
2835         effects[currentEffect]->SetBinContent(i, 0.01);
2836       }
2837       else if (i <= 81)
2838       {
2839         effects[currentEffect]->SetBinContent(i, 0.01 + 0.05 * (i - 61) / 20);
2840       }
2841     }
2842
2843 //     currentEffect++;
2844 //
2845 //     // material budget
2846 //     for (Int_t i=1; i<=81; ++i)
2847 //     {
2848 //       if (i < 5)
2849 //         effects[currentEffect]->SetBinContent(i, 0.05 - 0.01 * i);
2850 //       if (i > 51)
2851 //         effects[currentEffect]->SetBinContent(i, 0.05 * (i - 50) / 30);
2852 //     }
2853 //
2854     currentEffect++;
2855
2856   }
2857
2858   TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 500);
2859   canvas->SetRightMargin(0.05);
2860   canvas->SetTopMargin(0.05);
2861   //canvas->SetGridx();
2862   canvas->SetGridy();
2863   TLegend* legend = new TLegend(0.2, 0.4, 0.7, 0.4 + 0.5 * nEffects / 7);
2864   legend->SetFillColor(0);
2865   legend->SetTextSize(0.04);
2866   dummy->Draw();
2867   dummy->GetXaxis()->SetRangeUser(0, displayRange);
2868
2869   for (Int_t i=0; i<nEffects; ++i)
2870   {
2871     TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2872     /*current->Reset();
2873     for (Int_t j=0; j<nEffects-i; ++j)
2874       current->Add(effects[j]);*/
2875
2876     current->SetLineColor(colors[i]);
2877     current->SetLineStyle(styles[i]);
2878     current->SetLineWidth(widths[i]);
2879     //current->SetFillColor(colors[i]);
2880     current->SetMarkerColor(colors[i]);
2881     //current->SetMarkerStyle(markers[i]);
2882
2883     current->SetStats(kFALSE);
2884     current->GetYaxis()->SetRangeUser(0, 0.4);
2885     current->DrawCopy("SAME");
2886     legend->AddEntry(current, names[i]);
2887
2888     //TLatex* text = new TLatex(displayRange+2, current->GetBinContent(displayRange+1), names[i]);
2889     TLatex* text = new TLatex(displayRange+2, 0.1 - i * 0.02, names[i]);
2890     text->SetTextSize(0.04);
2891     text->SetTextColor(colors[i]);
2892     //text->Draw();
2893   }
2894
2895   // add total in square
2896   TH1* totalINEL = (TH1*) effects[0]->Clone("totalINEL");
2897   totalINEL->Reset();
2898   TH1* totalNSD = (TH1*) totalINEL->Clone("totalNSD");
2899
2900   for (Int_t i=0; i<nEffects; ++i)
2901   {
2902     //Printf("%d %f", i, effects[i]->GetBinContent(20));
2903     effects[i]->Multiply(effects[i]);
2904
2905     if (i != 2)
2906       totalINEL->Add(effects[i]);
2907     if (i != 1)
2908       totalNSD->Add(effects[i]);
2909   }
2910
2911   for (Int_t i=1; i<=totalINEL->GetNbinsX(); ++i)
2912   {
2913     totalINEL->SetBinError(i, 0);
2914     if (totalINEL->GetBinContent(i) > 0)
2915       totalINEL->SetBinContent(i, TMath::Min(sqrt(totalINEL->GetBinContent(i)), 1.0));
2916     totalNSD->SetBinError(i, 0);
2917     if (totalNSD->GetBinContent(i) > 0)
2918       totalNSD->SetBinContent(i, TMath::Min(sqrt(totalNSD->GetBinContent(i)), 1.0));
2919   }
2920
2921   //Printf("%f", total->GetBinContent(20));
2922
2923   totalINEL->SetMarkerStyle(5);
2924   totalINEL->SetMarkerColor(1);
2925   legend->AddEntry(totalINEL, "Total (INEL)", "P");
2926
2927   totalNSD->SetMarkerStyle(24);
2928   totalNSD->SetMarkerColor(2);
2929   legend->AddEntry(totalNSD, "Total (NSD)", "P");
2930
2931   Printf("total in bin 0 is INEL: %f NSD: %f", totalINEL->GetBinContent(1), totalNSD->GetBinContent(1));
2932   totalINEL->DrawCopy("SAME P"); //->SetBinContent(1, 0);
2933   totalNSD->DrawCopy("SAME P"); //->SetBinContent(1, 0);
2934
2935   legend->Draw();
2936
2937   canvas->SaveAs(canvas->GetName());
2938
2939   file = TFile::Open(Form("%s_syst_error.root", (tpc) ? "tpc" : "spd"), "RECREATE");
2940   totalINEL->Write("inel_1");
2941   totalNSD->Write("nsd_1");
2942   file->Close();
2943
2944   return (nsd) ? totalNSD : totalINEL;
2945 }
2946
2947 void finalPlot(Bool_t tpc = 0, Bool_t small = kFALSE)
2948 {
2949   loadlibs();
2950
2951   if (tpc)
2952     SetTPC();
2953
2954   //TH1* errorNSD = SystematicsSummary(tpc, 1);
2955
2956   TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", 800, 500);
2957   canvas->SetRightMargin(0.05);
2958   canvas->SetTopMargin(0.05);
2959   canvas->SetGridx();
2960   canvas->SetGridy();
2961
2962   legend = new TLegend(0.5, 0.6, 0.9, 0.8);
2963   legend->SetFillColor(0);
2964   legend->SetTextSize(0.04);
2965
2966   for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
2967   {
2968     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
2969     TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
2970     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2971
2972     DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
2973
2974     // normalize result
2975     //result->Scale(1.0 / result->Integral(2, displayRange));
2976
2977     result->GetXaxis()->SetRangeUser(0, displayRange);
2978     //result->SetBinContent(1, 0); result->SetBinError(1, 0);
2979     result->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (etaRange+1) * 0.5));
2980     result->SetMarkerStyle(0);
2981     result->SetLineColor(1);
2982     result->SetStats(kFALSE);
2983
2984     // systematic error
2985     TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
2986
2987     TH1* systError = (TH1*) result->Clone("systError");
2988     for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
2989       systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2990
2991     // change error drawing style
2992     systError->SetFillColor(15);
2993
2994     if (eventType == AliMultiplicityCorrection::kNSD)
2995     {
2996       result->SetLineColor(2);
2997       result->SetMarkerColor(2);
2998       result->SetMarkerStyle(5);
2999     }
3000
3001     canvas->cd();
3002     systError->DrawCopy(Form("E2 ][ %s", (eventType == AliMultiplicityCorrection::kINEL) ? "" : "SAME"));
3003     result->DrawCopy("SAME E ][");
3004     canvas->SetLogy();
3005
3006     legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic cross-section" : "NSD cross-section", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
3007   }
3008
3009   legend->Draw();
3010   /*
3011   //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
3012   TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
3013   text->SetFillColor(0);
3014   text->SetTextAlign(12);
3015   text->AddText("Systematic errors summed quadratically");
3016   text->AddText("0.6 million minimum bias events");
3017   text->AddText("corrected to inelastic events");
3018   text->Draw("B");
3019
3020   TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
3021   text2->SetFillColor(0);
3022   text2->SetTextAlign(12);
3023   text2->AddText("#sqrt{s} = 14 TeV");
3024   if (tpc)
3025   {
3026     text2->AddText("|#eta| < 0.9");
3027   }
3028   else
3029     text2->AddText("|#eta| < 2.0");
3030   text2->AddText("simulated data (PYTHIA)");
3031   text2->Draw("B");
3032
3033   if (tpc)
3034   {
3035     TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
3036     text3->SetNDC();
3037     text3->Draw();
3038   }
3039   else
3040   {
3041     TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
3042     text3->SetNDC();
3043     text3->Draw();
3044   }
3045
3046   // alice logo
3047   TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
3048   pad->Draw();
3049   pad->cd();
3050   TImage* img = TImage::Open("$HOME/alice.png");
3051   img->SetImageQuality(TAttImage::kImgBest);
3052   img->Draw();
3053
3054   canvas->Modified();
3055   */
3056
3057 /*  TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
3058   text->SetTextSize(0.04);
3059   text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
3060   text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
3061   text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
3062   text->Draw();*/
3063
3064
3065   canvas->SaveAs(canvas->GetName());
3066 }
3067
3068 void finalPlot2(Bool_t tpc = 0)
3069 {
3070   loadlibs();
3071
3072   if (tpc)
3073     SetTPC();
3074
3075   //TH1* errorNSD = SystematicsSummary(tpc, 1);
3076
3077   TCanvas* canvas = new TCanvas("finalPlot2.eps", "finalPlot2.eps", 800, 600);
3078   canvas->SetRightMargin(0.05);
3079   canvas->SetTopMargin(0.05);
3080   canvas->SetGridx();
3081   canvas->SetGridy();
3082   canvas->SetLogy();
3083
3084   legend = new TLegend(0.5, 0.7, 0.9, 0.9);
3085   legend->SetFillColor(0);
3086   legend->SetTextSize(0.04);
3087
3088   Int_t displayRanges[] = { 30, 45, 65 };
3089
3090   TH2* dummy = new TH2F("dummy", ";True multiplicity N_{ch};P(N_{ch})", 100, -0.5, displayRanges[2]+10, 1000, 5e-6, 5);
3091   dummy->SetStats(0);
3092   dummy->Draw();
3093
3094   TList list;
3095
3096   // get syst error
3097   file = TFile::Open(Form("%s_syst_error.root", (tpc) ? "tpc" : "spd"));
3098   TH1* totalINEL = (TH1*) file->Get("inel_1");
3099   TH1* totalNSD = (TH1*) file->Get("nsd_1");
3100
3101   Int_t count = 0;
3102   for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
3103   {
3104     AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
3105     //TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
3106     for (Int_t etaR = 0; etaR <= 2; etaR++)
3107     {
3108       TH1* result = mult->GetMultiplicityESDCorrected(etaR);
3109
3110       //DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
3111
3112       // normalize result
3113       result->Scale(TMath::Power(5, etaR) / result->Integral(1, displayRanges[etaR]));
3114
3115       result->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
3116       //result->SetBinContent(1, 0); result->SetBinError(1, 0);
3117       //result->SetTitle(Form(""));
3118       result->SetMarkerStyle(0);
3119       result->SetLineColor(1);
3120       result->SetLineWidth(2);
3121       //result->SetMarkerStyle(4);
3122       //result->SetStats(kFALSE);
3123
3124       // systematic error
3125       //TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
3126       TH1* error = (eventType == AliMultiplicityCorrection::kNSD) ? totalNSD : totalINEL;
3127
3128       TH1* systError = (TH1*) result->Clone("systError");
3129       // small hack until we have syst. errors for all eta ranges
3130       Float_t factor = 80.0 / displayRanges[etaR];
3131       for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
3132       {
3133         systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(factor * i));
3134       }
3135
3136       // change error drawing style
3137       systError->SetFillColor(15);
3138
3139       if (eventType == AliMultiplicityCorrection::kNSD)
3140       {
3141         result->SetLineColor(2);
3142         result->SetMarkerColor(2);
3143         result->SetMarkerStyle(5);
3144       }
3145
3146       canvas->cd();
3147       systError->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
3148       systError->DrawCopy("E2 ][ SAME");
3149       list.Add(result);
3150
3151       if (eventType == AliMultiplicityCorrection::kINEL)
3152       {
3153         TLatex* Tl = new TLatex;
3154         Tl->SetTextSize(0.04);
3155         //Tl->SetBit(TLatex::kTextNDC);
3156         Tl->SetTextAlign(32);
3157
3158         Float_t etaRangeArr[] = { 0.5, 1.0, 1.4 };
3159         TString tmpStr;
3160         tmpStr.Form("|#eta| < %.1f (x %d)", etaRangeArr[etaR], (Int_t) TMath::Power(5, etaR));
3161         if (etaR == 0)
3162           Tl->DrawLatex(15, result->GetBinContent(20), tmpStr);
3163         if (etaR == 1)
3164         {
3165           Tl->SetTextAlign(12);
3166           Tl->DrawLatex(40, result->GetBinContent(40), tmpStr);
3167         }
3168         if (etaR == 2)
3169         {
3170           Tl->SetTextAlign(12);
3171           Tl->DrawLatex(60, result->GetBinContent(50), tmpStr);
3172         }
3173       }
3174
3175       if (etaR == 0)
3176         legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic events" : "NSD events", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
3177
3178       count++;
3179     }
3180   }
3181
3182   for (Int_t i=0; i<list.GetEntries(); i++)
3183     ((TH1*) list.At(i))->DrawCopy("SAME E ][");
3184
3185   legend->Draw();
3186
3187   canvas->SaveAs(canvas->GetName());
3188 }
3189
3190 TMatrixD* NonInvertable()
3191 {
3192   const Int_t kSize = 5;
3193
3194   TMatrixD matrix(kSize, kSize);
3195   for (Int_t x=0; x<kSize; x++)
3196   {
3197     for (Int_t y=0; y<kSize; y++)
3198     {
3199       if (x == y)
3200       {
3201         if (x == 0 || x == kSize -1)
3202         {
3203           matrix(x, y) = 0.75;
3204         }
3205         else
3206           matrix(x, y) = 0.5;
3207       }
3208       else if (TMath::Abs(x - y) == 1)
3209       {
3210         matrix(x, y) = 0.25;
3211       }
3212     }
3213   }
3214
3215   matrix.Print();
3216
3217   //TMatrixD inverted(matrix);
3218   //inverted.Invert();
3219
3220   //inverted.Print();
3221
3222   return new TMatrixD(matrix);
3223 }
3224
3225 void BlobelUnfoldingExample()
3226 {
3227   const Int_t kSize = 20;
3228
3229   TMatrixD matrix(kSize, kSize);
3230   for (Int_t x=0; x<kSize; x++)
3231   {
3232     for (Int_t y=0; y<kSize; y++)
3233     {
3234       if (x == y)
3235       {
3236         if (x == 0 || x == kSize -1)
3237         {
3238           matrix(x, y) = 0.75;
3239         }
3240         else
3241           matrix(x, y) = 0.5;
3242       }
3243       else if (TMath::Abs(x - y) == 1)
3244       {
3245         matrix(x, y) = 0.25;
3246       }
3247     }
3248   }
3249
3250   //matrix.Print();
3251
3252   TMatrixD inverted(matrix);
3253   inverted.Invert();
3254
3255   //inverted.Print();
3256
3257   TH1F* inputDist = new TH1F("inputDist", ";t;Entries", kSize, -0.5, (Float_t) kSize - 0.5);
3258   TVectorD inputDistVector(kSize);
3259   TH1F* unfolded = inputDist->Clone("unfolded");
3260   TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
3261   measuredIdealDist->SetTitle(";m;Entries");
3262   TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
3263
3264   TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
3265   // norm: 1/(sqrt(2pi)sigma)
3266   gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
3267   //gaus->Print();
3268
3269   for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
3270   {
3271     Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
3272     inputDist->SetBinContent(x, value);
3273     inputDistVector(x-1) = value;
3274   }
3275
3276   TVectorD measuredDistIdealVector = matrix * inputDistVector;
3277
3278   for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
3279     measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
3280
3281   measuredDist->FillRandom(measuredIdealDist, 10000);
3282
3283   // fill error matrix before scaling
3284   TMatrixD covarianceMatrixMeasured(kSize, kSize);
3285   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3286     covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
3287
3288   TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
3289   //covarianceMatrix.Print();
3290
3291   TVectorD measuredDistVector(kSize);
3292   for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
3293     measuredDistVector(x-1) = measuredDist->GetBinContent(x);
3294
3295   TVectorD unfoldedVector = inverted * measuredDistVector;
3296   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
3297     unfolded->SetBinContent(x, unfoldedVector(x-1));
3298
3299   TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1200, 600);
3300   canvas->SetTopMargin(0.05);
3301   canvas->Divide(2, 1);
3302
3303   canvas->cd(1);
3304   canvas->cd(1)->SetLeftMargin(0.15);
3305   canvas->cd(1)->SetRightMargin(0.05);
3306   canvas->cd(1)->SetTopMargin(0.05);
3307   gPad->SetGridx();
3308   gPad->SetGridy();
3309   measuredDist->GetYaxis()->SetRangeUser(-600, 2799);
3310   measuredDist->GetYaxis()->SetTitleOffset(1.9);
3311   measuredDist->SetStats(0);
3312   measuredDist->DrawCopy();
3313   gaus->Draw("SAME");
3314
3315   canvas->cd(2);
3316   canvas->cd(2)->SetLeftMargin(0.15);
3317   canvas->cd(2)->SetRightMargin(0.05);
3318   canvas->cd(2)->SetTopMargin(0.05);
3319   gPad->SetGridx();
3320   gPad->SetGridy();
3321   unfolded->GetYaxis()->SetRangeUser(-600, 2799);
3322   unfolded->GetYaxis()->SetTitleOffset(1.9);
3323   unfolded->SetStats(0);
3324   unfolded->DrawCopy();
3325   gaus->Draw("SAME");
3326
3327   canvas->SaveAs("BlobelUnfoldingExample.eps");
3328
3329   return;
3330
3331   // now unfold this with Bayesian
3332   loadlibs();
3333
3334   // fill a multiplicity object
3335   mult = new AliMultiplicityCorrection("mult", "mult");
3336   for (Int_t x=0; x<kSize; x++)
3337   {
3338     mult->GetMultiplicityVtx(0)->SetBinContent(1, x+1, inputDistVector(x));
3339     mult->GetMultiplicityESD(0)->SetBinContent(1, x+1, measuredDistVector(x)*10000);
3340     for (Int_t y=0; y<kSize; y++)
3341       mult->GetCorrelation(0)->SetBinContent(1, x+1, y+1, matrix(x, y));
3342   }
3343
3344   //mult->DrawHistograms();
3345
3346   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 0);
3347   //mult->SetCreateBigBin(kFALSE);
3348   mult->ApplyMinuitFit(0, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); //hist2->ProjectionY("mymchist"));
3349
3350   //mult->ApplyBayesianMethod(0, kFALSE, AliMultiplicityCorrection::kTrVtx, 0, -1, 0, kFALSE);
3351
3352   mult->DrawComparison("BlobelExample", 0, kFALSE, kTRUE, mult->GetMultiplicityVtx(0)->ProjectionY("mcmchist", 1, mult->GetMultiplicityVtx(0)->GetNbinsX()));
3353 }
3354
3355 void TestWithDiagonalMatrix()
3356 {
3357   const Int_t kSize = 20;
3358
3359   TMatrixD matrix(kSize, kSize);
3360   for (Int_t x=0; x<kSize; x++)
3361     matrix(x, x) = 1;
3362
3363   if (1)
3364   {
3365   for (Int_t x=0; x<kSize; x++)
3366   {
3367     for (Int_t y=0; y<kSize; y++)
3368     {
3369       if (x == y)
3370       {
3371         if (x == 0 || x == kSize -1)
3372         {
3373           matrix(x, y) = 0.75;
3374         }
3375         else
3376           matrix(x, y) = 0.5;
3377       }
3378       else if (TMath::Abs(x - y) == 1)
3379       {
3380         matrix(x, y) = 0.25;
3381       }
3382     }
3383   }
3384   }
3385
3386   //matrix.Print();
3387
3388   TH1F* inputDist = new TH1F("inputDist", ";t;Entries", kSize, -0.5, (Float_t) kSize - 0.5);
3389   TVectorD inputDistVector(kSize);
3390
3391   TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
3392   measuredIdealDist->SetTitle(";m;Entries");
3393   TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
3394
3395   TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
3396   // norm: 1/(sqrt(2pi)sigma)
3397   gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
3398   //gaus->Print();
3399
3400   for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
3401   {
3402     Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
3403     inputDist->SetBinContent(x, value);
3404     inputDistVector(x-1) = value;
3405   }
3406
3407   TVectorD measuredDistIdealVector = matrix * inputDistVector;
3408
3409   for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
3410     measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
3411
3412   // randomize
3413   //measuredDist->FillRandom(measuredIdealDist, 10000);
3414   measuredDist = measuredIdealDist;
3415   measuredDist->Sumw2();
3416
3417   for (Int_t x=0; x<kSize; x++)
3418     Printf("bin %d %.2f +- %.2f", x+1, measuredDist->GetBinContent(x+1), measuredDist->GetBinError(x+1));
3419
3420   // now unfold this
3421   loadlibs();
3422
3423   // fill a multiplicity object
3424   mult = new AliMultiplicityCorrection("mult", "mult");
3425   for (Int_t x=0; x<kSize; x++)
3426   {
3427     mult->GetMultiplicityVtx(0)->SetBinContent(1, x+1, inputDistVector(x));
3428     mult->GetMultiplicityESD(0)->SetBinContent(1, x+1, measuredDist->GetBinContent(x+1));
3429     for (Int_t y=0; y<kSize; y++)
3430       mult->GetCorrelation(0)->SetBinContent(1, x+1, y+1, matrix(x, y));
3431   }
3432
3433   //mult->DrawHistograms();
3434
3435   AliUnfolding::SetNbins(20, 20);
3436   AliUnfolding::SetChi2Regularization(AliUnfolding::kNone, 0);
3437   AliUnfolding::SetChi2Regularization(AliUnfolding::kPol1, 1e-14);
3438   //mult->SetCreateBigBin(kFALSE);
3439   mult->ApplyMinuitFit(0, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); //hist2->ProjectionY("mymchist"));
3440
3441   //mult->ApplyBayesianMethod(0, kFALSE, AliMultiplicityCorrection::kTrVtx, 0, -1, 0, kFALSE);
3442
3443   mult->DrawComparison("TestWithDiagonalMatrix", 0, kFALSE, kTRUE, mult->GetMultiplicityVtx(0)->ProjectionY("mcmchist", 1, mult->GetMultiplicityVtx(0)->GetNbinsX()));
3444 }
3445
3446
3447 void E735Fit()
3448 {
3449   TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
3450   fCurrentESD->Sumw2();
3451
3452   // Open the input stream
3453   ifstream in;
3454   in.open("e735data.txt");
3455
3456   while(in.good())
3457   {
3458     Float_t x, y, ye;
3459     in >> x >> y >> ye;
3460
3461     //Printf("%f %f %f", x, y, ye);
3462     fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
3463     fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
3464   }
3465
3466   in.close();
3467
3468   //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
3469
3470   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
3471
3472   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])");
3473   func->SetParNames("scaling", "averagen", "k");
3474   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
3475   func->SetParLimits(1, 0.001, 1000);
3476   func->SetParLimits(2, 0.001, 1000);
3477   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
3478
3479   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
3480   lognormal->SetParNames("scaling", "mean", "sigma");
3481   lognormal->SetParameters(1, 1, 1);
3482   lognormal->SetParLimits(0, 0, 10);
3483   lognormal->SetParLimits(1, 0, 100);
3484   lognormal->SetParLimits(2, 1e-3, 10);
3485
3486   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
3487   fCurrentESD->SetStats(kFALSE);
3488   fCurrentESD->SetMarkerStyle(0);
3489   fCurrentESD->SetLineColor(1);
3490   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
3491   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
3492   fCurrentESD->Draw("");
3493   fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
3494   fCurrentESD->Fit(func, "0", "", 0, 150);
3495   func->SetRange(0, 250);
3496   func->Draw("SAME");
3497   printf("chi2 = %f\n", func->GetChisquare());
3498
3499   fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
3500   lognormal->SetLineColor(2);
3501   lognormal->SetLineStyle(2);
3502   lognormal->SetRange(0, 250);
3503   lognormal->Draw("SAME");
3504
3505   gPad->SetLogy();
3506
3507   canvas->SaveAs("E735Fit.eps");
3508 }
3509
3510 void DifferentSamples()
3511 {
3512   loadlibs();
3513
3514   Int_t n = 2;
3515   const char* filesChi2[] = { "chi2_100k_1.root", "chi2_100k_2.root" };
3516   const char* filesBayesian[] = { "bayesian_100k_1.root", "bayesian_100k_2.root" };
3517
3518   TCanvas* canvas = new TCanvas("DifferentSamples", "DifferentSamples", 1200, 600);
3519   canvas->Divide(2, 1);
3520
3521   legend = new TLegend(0.15, 0.7, 0.65, 0.9);
3522   legend->SetFillColor(0);
3523   legend->SetTextSize(0.04);
3524
3525   for (Int_t i=0; i<n; i++)
3526   {
3527     AliMultiplicityCorrection* chi2 = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3528     TFile::Open(filesChi2[i]);
3529     chi2->LoadHistograms("Multiplicity");
3530
3531     AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3532     TFile::Open(filesBayesian[i]);
3533     bayesian->LoadHistograms("Multiplicity");
3534
3535     chi2Hist = chi2->GetMultiplicityESDCorrected(etaRange);
3536     bayesianHist = bayesian->GetMultiplicityESDCorrected(etaRange);
3537
3538     mc = chi2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
3539
3540     // normalize and divide
3541     chi2Hist->Scale(1.0 / chi2Hist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3542     bayesianHist->Scale(1.0 / bayesianHist->Integral(2, displayRange+1) * mc->Integral(2, displayRange));
3543
3544     chi2Hist->Divide(mc, chi2Hist);
3545     bayesianHist->Divide(mc, bayesianHist);
3546
3547     canvas->cd(i+1);
3548     gPad->SetTopMargin(0.05);
3549     gPad->SetRightMargin(0.05);
3550     //gPad->SetLeftMargin(0.12);
3551     gPad->SetGridx();
3552     gPad->SetGridy();
3553
3554     chi2Hist->GetXaxis()->SetRangeUser(0, displayRange);
3555     chi2Hist->GetYaxis()->SetTitleOffset(1.3);
3556     chi2Hist->SetStats(0);
3557     chi2Hist->SetTitle(Form(";%s;MC / unfolded", GetMultLabel()));
3558     chi2Hist->GetYaxis()->SetRangeUser(0.2, 1.8);
3559     chi2Hist->Draw("HIST");
3560
3561     for (Int_t x=1; x<=bayesianHist->GetNbinsX(); x++)
3562       bayesianHist->SetBinError(x, 1e-6);
3563
3564     bayesianHist->SetLineColor(2);
3565     bayesianHist->SetMarkerColor(2);
3566     bayesianHist->SetMarkerStyle(5);
3567     bayesianHist->Draw("HIST E SAME");
3568
3569     if (i == 0)
3570     {
3571       legend->AddEntry(chi2Hist, "#chi^{2}-minimization", "L");
3572       legend->AddEntry(bayesianHist, "Bayesian unfolding", "LP");
3573     }
3574     legend->Draw();
3575   }
3576
3577   canvas->SaveAs("DifferentSamples.eps");
3578 }
3579
3580 void PileUp()
3581 {
3582   loadlibs();
3583
3584   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
3585   hist2d = mult->GetMultiplicityMC(etaRange, AliMultiplicityCorrection::kINEL);
3586   mult1 = hist2d->ProjectionY("mult1", 1, hist2d->GetNbinsX());
3587
3588   conv = (TH1*) mult1->Clone("conv");
3589   conv->Reset();
3590
3591   mult1->Scale(1.0 / mult1->Integral());
3592
3593   for (Int_t i=1; i<=mult1->GetNbinsX(); i++)
3594     for (Int_t j=1; j<=mult1->GetNbinsX(); j++)
3595       conv->Fill(mult1->GetBinCenter(i)+mult1->GetBinCenter(j), mult1->GetBinContent(i) * mult1->GetBinContent(j));
3596
3597   conv->Scale(1.0 / conv->Integral());
3598
3599   c = new TCanvas("c", "c", 800, 500);
3600   gPad->SetLogy();
3601   gPad->SetTopMargin(0.05);
3602   gPad->SetRightMargin(0.05);
3603   mult1->SetTitle(Form(";%s;Probability", GetMultLabel()));
3604   mult1->SetStats(0);
3605   gPad->SetGridx();
3606   gPad->SetGridy();
3607   mult1->Draw();
3608   mult1->GetYaxis()->SetRangeUser(1e-7, 2 * mult1->GetMaximum());
3609   mult1->GetXaxis()->SetRangeUser(0, displayRange);
3610   mult1->GetXaxis()->SetTitleOffset(1.15);
3611   conv->SetLineColor(2);
3612   conv->SetMarkerColor(2);
3613   conv->SetMarkerStyle(5);
3614   conv->DrawCopy("SAME P");
3615
3616   conv->Scale(0.00058);
3617   conv->DrawCopy("SAME P");
3618
3619   legend = new TLegend(0.73, 0.73, 0.93, 0.93);
3620   legend->SetFillColor(0);
3621   legend->SetTextSize(0.04);
3622   legend->AddEntry(mult1, "1 collision");
3623   legend->AddEntry(conv, "2 collisions", "P");
3624   legend->Draw();
3625
3626   c->SaveAs("pileup.eps");
3627
3628   new TCanvas;
3629   conv->Divide(mult1);
3630   conv->Draw();
3631 }
3632
3633 void TestErrorDetermination(Int_t nRandomizations)
3634 {
3635   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);
3636   func->SetParNames("scaling", "averagen", "k");
3637   func->SetParameters(1, 15, 2);
3638
3639   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);
3640   func2->SetParNames("scaling", "averagen", "k");
3641   func2->SetParLimits(0, 0.5, 2);
3642   func2->SetParLimits(1, 1, 50);
3643   func2->SetParLimits(2, 1, 10);
3644   func2->SetParameters(1, 15, 2);
3645   //func2->FixParameter(0, 1);
3646
3647   //new TCanvas; func->Draw("L");
3648
3649   hist1 = new TH1F("hist1", "", 100, 0.5, 100.5);
3650   hist2 = new TH1F("hist2", "", 100, 0.5, 100.5);
3651   hist1->Sumw2();
3652
3653   TH1* params[3];
3654   params[0] = new TH1F("param_0", Form("param_%d", 0), 100, 0.95, 1.05);
3655   params[1] = new TH1F("param_1", Form("param_%d", 1), 100, 14, 16);
3656   params[2] = new TH1F("param_2", Form("param_%d", 2), 100, 1.8, 2.2);
3657
3658   const Int_t nTries = 1000;
3659   for (Int_t i=0; i<nTries; i++)
3660   {
3661     hist1->Reset();
3662
3663     if (nRandomizations == 1)
3664     {
3665       hist1->FillRandom("nbd", 10000);
3666     }
3667     else if (nRandomizations == 2)
3668     {
3669       hist2->Reset();
3670       hist2->FillRandom("nbd", 10000);
3671       hist1->FillRandom(hist2, 10000);
3672     }
3673     else if (nRandomizations == 3)
3674     {
3675       hist2->Reset();
3676       hist1->FillRandom("nbd", 10000);
3677       hist2->FillRandom(hist1, 10000);
3678       hist1->Reset();
3679       hist1->FillRandom(hist2, 10000);
3680     }
3681     else
3682       return;
3683
3684     //new TCanvas; hist1->Draw();
3685
3686     hist1->Scale(1.0 / hist1->Integral());
3687     hist1->Fit(func2, "NQ");
3688     hist1->Fit(func2, "NQ");
3689     for (Int_t j=0; j<3; j++)
3690       params[j]->Fill(func2->GetParameter(j));
3691   }
3692
3693   for (Int_t j=0; j<3; j++)
3694   {
3695     new TCanvas; params[j]->Draw();
3696     params[j]->Fit("gaus");
3697     Printf("sigma of param %d if %f", j, ((TF1*) params[j]->FindObject("gaus"))->GetParameter(2));
3698   }
3699 }
3700
3701 void DrawRawDistributions(const char* fileName = "multiplicityESD.root")
3702 {
3703   loadlibs();
3704
3705   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileName);
3706
3707   c = new TCanvas("c", "c", 600, 600);
3708
3709   dummy = new TH2F("dummy", ";measured multiplicity", 100, -0.5, 149.5, 100, 0.5, 4e4);
3710   dummy->SetStats(0);
3711   dummy->Draw();
3712   gPad->SetGridx();
3713   gPad->SetGridy();
3714   gPad->SetLogy();
3715
3716   Int_t colors[] = { 1, 2, 4 };
3717
3718   for (Int_t i=2; i>=0; i--)
3719   {
3720     hist = mult->GetMultiplicityESD(i)->ProjectionY();
3721
3722     hist->SetLineColor(colors[i]);
3723     hist->DrawCopy("SAME");
3724   }
3725
3726
3727 }
3728
3729 void FindUnfoldedLimit()
3730 {
3731   loadlibs();
3732
3733   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
3734   AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open("multiplicityESD.root");
3735
3736   TH1* hist = mult->GetCorrelation(etaRange);
3737   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
3738   {
3739     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
3740     {
3741       hist->SetBinContent(0, y, z, 0);
3742       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
3743     }
3744   }
3745   TH2* corr = (TH2*) ((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy");
3746
3747   TH1* esd_proj = esd->GetMultiplicityESD(etaRange)->ProjectionY("esd_proj", 1, esd->GetMultiplicityESD(etaRange)->GetNbinsX());
3748   //esd_proj = corr->ProjectionY("esd_proj", 1, corr->GetNbinsX());
3749
3750   new TCanvas; corr->Draw("COLZ");
3751   new TCanvas; esd_proj->DrawCopy();
3752
3753   TH1* percentage = (TH1*) (esd_proj->Clone("percentage"));
3754   percentage->Reset();
3755
3756   for (Int_t i=1; i<=esd_proj->GetNbinsX(); i++)
3757     if (esd_proj->GetBinContent(i) > 0)
3758       esd_proj->SetBinContent(i, 1);
3759
3760   for (Int_t i=1; i<=percentage->GetNbinsX(); i++)
3761   {
3762     TH1* binResponse = corr->ProjectionY("proj", i, i);
3763     if (binResponse->Integral() <= 0)
3764       continue;
3765     binResponse->Scale(1.0 / binResponse->Integral());
3766     binResponse->Multiply(esd_proj);
3767     //new TCanvas; binResponse->Draw();
3768     percentage->SetBinContent(i, binResponse->Integral());
3769     //return;
3770   }
3771
3772   new TCanvas; percentage->Draw();
3773   new TCanvas;
3774   mc = esd->GetMultiplicityVtx(etaRange)->ProjectionY("mc", 1, esd->GetMultiplicityVtx(etaRange)->GetNbinsX());
3775   mc->SetLineColor(2);
3776   mc->Draw("");
3777 }
3778
3779 void CompareUnfoldedWithMC()
3780 {
3781   loadlibs();
3782
3783   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("unfolded.root");
3784
3785   //mult->GetMultiplicityESDCorrected(etaRange)->Scale(1.0 / mult->GetMultiplicityESDCorrected(etaRange)->Integral());
3786   mult->GetMultiplicityESDCorrected(etaRange)->SetTitle(";multiplicity;events");
3787   mult->GetMultiplicityESDCorrected(etaRange)->SetStats(0);
3788   mult->GetMultiplicityESDCorrected(etaRange)->Draw();
3789   //mcHist->Scale(1.0 / mcHist->Integral());
3790   mcHist = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
3791   mcHist->SetLineColor(2);
3792   mcHist->Draw("SAME");
3793   gPad->SetLogy();
3794 }