PWG0/dNdEta/plotsMultiplicity.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_2M.root";
  34 const char* measuredFile   = "multiplicityMC_1M_3.root";
  35 Int_t etaRange = 3;
  36 Int_t displayRange = 200; // axis range
  37 Int_t ratioRange = 151;   // range to calculate difference
  38 Int_t longDisplayRange = 200;
  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 SetTPC()
  45 {
  46   correctionFile = correctionFileTPC;
  47   measuredFile = measuredFileTPC;
  48   etaRange = etaRangeTPC;
  49   displayRange = 100;
  50   ratioRange = 76;
  51   longDisplayRange = 100;
  52 }
  53
  54 void Smooth(TH1* hist, Int_t windowWidth = 20)
  55 {
  56   TH1* clone = (TH1*) hist->Clone("clone");
  57   for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
  58   {
  59     Int_t min = TMath::Max(2, bin-windowWidth);
  60     Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
  61     Float_t average = clone->Integral(min, max) / (max - min + 1);
  62
  63     hist->SetBinContent(bin, average);
  64     hist->SetBinError(bin, 0);
  65   }
  66
  67   delete clone;
  68 }
  69
  70 void responseMatrixPlot()
  71 {
  72   gSystem->Load("libPWG0base");
  73
  74   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
  75
  76   TFile::Open(correctionFile);
  77   mult->LoadHistograms("Multiplicity");
  78
  79   TH1* hist = mult->GetCorrelation(etaRange)->Project3D("zy");
  80   hist->SetStats(kFALSE);
  81
  82   hist->SetTitle(";true multiplicity;measured multiplicity;Entries");
  83   hist->GetXaxis()->SetRangeUser(0, longDisplayRange);
  84   hist->GetYaxis()->SetRangeUser(0, longDisplayRange);
  85
  86   TCanvas* canvas = new TCanvas("c1", "c1", 800, 600);
  87   canvas->SetRightMargin(0.15);
  88   canvas->SetTopMargin(0.05);
  89
  90   gPad->SetLogz();
  91   hist->Draw("COLZ");
  92
  93   canvas->SaveAs("responsematrix.eps");
  94 }
  95
  96 TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
  97 {
  98   // normalize unfolded result to mc hist
  99   result->Scale(1.0 / result->Integral(2, 200));
 100   result->Scale(mcHist->Integral(2, 200));
 101
 102   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
 103   canvas->Range(0, 0, 1, 1);
 104
 105   TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
 106   pad1->Draw();
 107
 108   TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
 109   pad2->Draw();
 110
 111   pad1->SetRightMargin(0.05);
 112   pad2->SetRightMargin(0.05);
 113
 114   // no border between them
 115   pad1->SetBottomMargin(0);
 116   pad2->SetTopMargin(0);
 117
 118   pad1->cd();
 119
 120   mcHist->GetXaxis()->SetLabelSize(0.06);
 121   mcHist->GetYaxis()->SetLabelSize(0.06);
 122   mcHist->GetXaxis()->SetTitleSize(0.06);
 123   mcHist->GetYaxis()->SetTitleSize(0.06);
 124   mcHist->GetYaxis()->SetTitleOffset(0.6);
 125
 126   mcHist->GetXaxis()->SetRangeUser(0, displayRange);
 127
 128   mcHist->SetTitle(";true multiplicity;Entries");
 129   mcHist->SetStats(kFALSE);
 130
 131   mcHist->DrawCopy("HIST E");
 132   gPad->SetLogy();
 133
 134   result->SetLineColor(2);
 135   result->DrawCopy("SAME HISTE");
 136
 137   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
 138   legend->AddEntry(mcHist, "true distribution");
 139   legend->AddEntry(result, "unfolded distribution");
 140   legend->SetFillColor(0);
 141   legend->Draw();
 142
 143   pad2->cd();
 144   pad2->SetBottomMargin(0.15);
 145
 146   // calculate ratio
 147   mcHist->Sumw2();
 148   TH1* ratio = (TH1*) mcHist->Clone("ratio");
 149   result->Sumw2();
 150   ratio->Divide(ratio, result, 1, 1, "");
 151   ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
 152   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
 153
 154   ratio->DrawCopy();
 155
 156   // get average of ratio
 157   Float_t sum = 0;
 158   for (Int_t i=2; i<=ratioRange; ++i)
 159   {
 160     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
 161   }
 162   sum /= ratioRange-1;
 163
 164   printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
 165
 166   TLine* line = new TLine(0, 1, displayRange, 1);
 167   line->SetLineWidth(2);
 168   line->Draw();
 169
 170   line = new TLine(0, 1.1, displayRange, 1.1);
 171   line->SetLineWidth(2);
 172   line->SetLineStyle(2);
 173   line->Draw();
 174   line = new TLine(0, 0.9, displayRange, 0.9);
 175   line->SetLineWidth(2);
 176   line->SetLineStyle(2);
 177   line->Draw();
 178
 179   canvas->Modified();
 180
 181   canvas->SaveAs(epsName);
 182
 183   return canvas;
 184 }
 185
 186 TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
 187 {
 188   // draws the 3 plots in the upper plot
 189   // draws the ratio between result1 and result2 in the lower plot
 190
 191   // normalize unfolded result to mc hist
 192   result1->Scale(1.0 / result1->Integral(2, 200));
 193   result1->Scale(mcHist->Integral(2, 200));
 194   result2->Scale(1.0 / result2->Integral(2, 200));
 195   result2->Scale(mcHist->Integral(2, 200));
 196
 197   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
 198   canvas->Range(0, 0, 1, 1);
 199
 200   TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
 201   pad1->Draw();
 202
 203   TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
 204   pad2->Draw();
 205
 206   pad1->SetRightMargin(0.05);
 207   pad2->SetRightMargin(0.05);
 208
 209   // no border between them
 210   pad1->SetBottomMargin(0);
 211   pad2->SetTopMargin(0);
 212
 213   pad1->cd();
 214
 215   mcHist->GetXaxis()->SetLabelSize(0.06);
 216   mcHist->GetYaxis()->SetLabelSize(0.06);
 217   mcHist->GetXaxis()->SetTitleSize(0.06);
 218   mcHist->GetYaxis()->SetTitleSize(0.06);
 219   mcHist->GetYaxis()->SetTitleOffset(0.6);
 220
 221   mcHist->GetXaxis()->SetRangeUser(0, displayRange);
 222
 223   mcHist->SetTitle(";true multiplicity;Entries");
 224   mcHist->SetStats(kFALSE);
 225
 226   mcHist->DrawCopy("HIST E");
 227   gPad->SetLogy();
 228
 229   result1->SetLineColor(2);
 230   result1->DrawCopy("SAME HISTE");
 231
 232   result2->SetLineColor(4);
 233   result2->DrawCopy("SAME HISTE");
 234
 235   TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
 236   legend->AddEntry(mcHist, "true distribution");
 237   legend->AddEntry(result1, "unfolded distribution (syst)");
 238   legend->AddEntry(result2, "unfolded distribution (normal)");
 239   legend->SetFillColor(0);
 240   legend->Draw();
 241
 242   pad2->cd();
 243   pad2->SetBottomMargin(0.15);
 244
 245   result1->GetXaxis()->SetLabelSize(0.06);
 246   result1->GetYaxis()->SetLabelSize(0.06);
 247   result1->GetXaxis()->SetTitleSize(0.06);
 248   result1->GetYaxis()->SetTitleSize(0.06);
 249   result1->GetYaxis()->SetTitleOffset(0.6);
 250
 251   result1->GetXaxis()->SetRangeUser(0, displayRange);
 252
 253   result1->SetTitle(";true multiplicity;Entries");
 254   result1->SetStats(kFALSE);
 255
 256   // calculate ratio
 257   result1->Sumw2();
 258   TH1* ratio = (TH1*) result1->Clone("ratio");
 259   result2->Sumw2();
 260   ratio->Divide(ratio, result2, 1, 1, "");
 261   ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
 262   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
 263
 264   ratio->DrawCopy();
 265
 266   // get average of ratio
 267   Float_t sum = 0;
 268   for (Int_t i=2; i<=ratioRange; ++i)
 269   {
 270     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
 271   }
 272   sum /= ratioRange-1;
 273
 274   printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
 275
 276   TLine* line = new TLine(0, 1, displayRange, 1);
 277   line->SetLineWidth(2);
 278   line->Draw();
 279
 280   line = new TLine(0, 1.1, displayRange, 1.1);
 281   line->SetLineWidth(2);
 282   line->SetLineStyle(2);
 283   line->Draw();
 284   line = new TLine(0, 0.9, displayRange, 0.9);
 285   line->SetLineWidth(2);
 286   line->SetLineStyle(2);
 287   line->Draw();
 288
 289   canvas->Modified();
 290
 291   canvas->SaveAs(epsName);
 292
 293   return canvas;
 294 }
 295
 296 TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE)
 297 {
 298   // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
 299   // a systematic effect
 300
 301   // normalize results
 302   result->Scale(1.0 / result->Integral(2, 200));
 303
 304   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
 305   canvas->SetTopMargin(0.05);
 306   canvas->SetRightMargin(0.05);
 307
 308   result->GetXaxis()->SetRangeUser(0, displayRange);
 309   result->GetYaxis()->SetRangeUser(0.55, 1.45);
 310   result->SetStats(kFALSE);
 311
 312   // to get the axis how we want it
 313   TH1* dummy = (TH1*) result->Clone("dummy");
 314   dummy->Reset();
 315   dummy->SetTitle(";true multiplicity;Ratio");
 316   dummy->DrawCopy();
 317   delete dummy;
 318
 319   Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
 320
 321   TLegend* legend = new TLegend(0.2, 0.75, 0.35, 0.95);
 322   legend->SetFillColor(0);
 323
 324   for (Int_t n=0; n<nResultSyst; ++n)
 325   {
 326     resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
 327
 328     // calculate ratio
 329     TH1* ratio = (TH1*) result->Clone("ratio");
 330     ratio->Divide(ratio, resultSyst[n], 1, 1, "");
 331     ratio->GetXaxis()->SetRangeUser(1, displayRange);
 332
 333     if (firstMarker)
 334       ratio->SetMarkerStyle(5);
 335
 336     ratio->SetLineColor(colors[n / 2]);
 337     if ((n % 2))
 338       ratio->SetLineStyle(2);
 339
 340     TString drawStr("SAME HIST");
 341     if (n == 0 && firstMarker)
 342       drawStr = "SAME P";
 343     if (errors)
 344       drawStr += " E";
 345
 346     ratio->DrawCopy(drawStr);
 347
 348     if (legendStrings && legendStrings[n])
 349       legend->AddEntry(ratio, legendStrings[n]);
 350
 351     // get average of ratio
 352     Float_t sum = 0;
 353     for (Int_t i=2; i<=ratioRange; ++i)
 354       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
 355     sum /= ratioRange-1;
 356
 357     printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
 358   }
 359
 360   if (legendStrings)
 361     legend->Draw();
 362
 363   TLine* line = new TLine(0, 1, displayRange, 1);
 364   line->SetLineWidth(2);
 365   line->Draw();
 366
 367   line = new TLine(0, 1.1, displayRange, 1.1);
 368   line->SetLineWidth(2);
 369   line->SetLineStyle(2);
 370   line->Draw();
 371   line = new TLine(0, 0.9, displayRange, 0.9);
 372   line->SetLineWidth(2);
 373   line->SetLineStyle(2);
 374   line->Draw();
 375
 376   canvas->SaveAs(epsName);
 377   canvas->SaveAs(Form("%s.gif", epsName.Data()));
 378
 379   return canvas;
 380 }
 381
 382 TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
 383 {
 384   // draws the ratios of each mc to the corresponding result
 385
 386   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
 387   canvas->SetRightMargin(0.05);
 388   canvas->SetTopMargin(0.05);
 389
 390   for (Int_t n=0; n<nResultSyst; ++n)
 391   {
 392     // normalize
 393     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
 394     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
 395
 396     result[n]->GetXaxis()->SetRangeUser(0, displayRange);
 397     result[n]->SetStats(kFALSE);
 398
 399     // calculate ratio
 400     TH1* ratio = (TH1*) result[n]->Clone("ratio");
 401     ratio->Divide(mc[n], ratio, 1, 1, "B");
 402
 403     // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
 404     ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
 405
 406     if (smooth)
 407       Smooth(ratio);
 408
 409     ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
 410     ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
 411
 412     if (dashed)
 413     {
 414       ratio->SetLineColor((n/2)+1);
 415       ratio->SetLineStyle((n%2)+1);
 416     }
 417     else
 418       ratio->SetLineColor(n+1);
 419
 420     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
 421
 422     // get average of ratio
 423     Float_t sum = 0;
 424     for (Int_t i=2; i<=ratioRange; ++i)
 425       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
 426     sum /= ratioRange-1;
 427
 428     printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
 429   }
 430
 431   TLine* line = new TLine(0, 1, displayRange, 1);
 432   line->SetLineWidth(2);
 433   line->Draw();
 434
 435   line = new TLine(0, 1.1, displayRange, 1.1);
 436   line->SetLineWidth(2);
 437   line->SetLineStyle(2);
 438   line->Draw();
 439   line = new TLine(0, 0.9, displayRange, 0.9);
 440   line->SetLineWidth(2);
 441   line->SetLineStyle(2);
 442   line->Draw();
 443
 444   canvas->Modified();
 445
 446   canvas->SaveAs(epsName);
 447   canvas->SaveAs(Form("%s.gif", epsName.Data()));
 448
 449   return canvas;
 450 }
 451
 452 TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
 453 {
 454   // draws the ratios of each mc to the corresponding result
 455   // deducts from each ratio the ratio of mcBase / resultBase
 456
 457   // normalize
 458   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
 459   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
 460
 461   // calculate ratio
 462   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
 463   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
 464
 465   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
 466   canvas->SetRightMargin(0.05);
 467   canvas->SetTopMargin(0.05);
 468
 469   for (Int_t n=0; n<nResultSyst; ++n)
 470   {
 471     // normalize
 472     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
 473     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
 474
 475     result[n]->GetXaxis()->SetRangeUser(0, displayRange);
 476     result[n]->SetStats(kFALSE);
 477
 478     // calculate ratio
 479     TH1* ratio = (TH1*) result[n]->Clone("ratio");
 480     ratio->Divide(mc[n], ratio, 1, 1, "B");
 481     ratio->Add(ratioBase, -1);
 482
 483     ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
 484     ratio->GetYaxis()->SetRangeUser(-1, 1);
 485     ratio->SetLineColor(n+1);
 486     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
 487
 488     // get average of ratio
 489     Float_t sum = 0;
 490     for (Int_t i=2; i<=ratioRange; ++i)
 491       sum += TMath::Abs(ratio->GetBinContent(i));
 492     sum /= ratioRange-1;
 493
 494     printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, ratioRange, sum);
 495   }
 496
 497   TLine* line = new TLine(0, 0, displayRange, 0);
 498   line->SetLineWidth(2);
 499   line->Draw();
 500
 501   line = new TLine(0, 0.1, displayRange, 0.1);
 502   line->SetLineWidth(2);
 503   line->SetLineStyle(2);
 504   line->Draw();
 505   line = new TLine(0, -0.1, displayRange, -0.1);
 506   line->SetLineWidth(2);
 507   line->SetLineStyle(2);
 508   line->Draw();
 509
 510   canvas->Modified();
 511
 512   canvas->SaveAs(epsName);
 513   canvas->SaveAs(Form("%s.gif", epsName.Data()));
 514
 515   return canvas;
 516 }
 517
 518 TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
 519 {
 520   // draws the ratios of each mc to the corresponding result
 521   // deducts from each ratio the ratio of mcBase / resultBase
 522   // smoothens the ratios by a sliding window
 523
 524   // normalize
 525   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
 526   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
 527
 528   // calculate ratio
 529   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
 530   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
 531
 532   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
 533   canvas->SetRightMargin(0.05);
 534   canvas->SetTopMargin(0.05);
 535
 536   for (Int_t n=0; n<nResultSyst; ++n)
 537   {
 538     // normalize
 539     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
 540     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
 541
 542     result[n]->GetXaxis()->SetRangeUser(0, displayRange);
 543     result[n]->SetStats(kFALSE);
 544
 545     // calculate ratio
 546     TH1* ratio = (TH1*) result[n]->Clone("ratio");
 547     ratio->Divide(mc[n], ratio, 1, 1, "B");
 548     ratio->Add(ratioBase, -1);
 549
 550     //new TCanvas; ratio->DrawCopy();
 551     // clear 0 bin
 552     ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
 553
 554     Smooth(ratio);
 555
 556     //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
 557
 558     canvas->cd();
 559     ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
 560     ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
 561     ratio->SetLineColor((n / 2)+1);
 562     ratio->SetLineStyle((n % 2)+1);
 563     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
 564
 565     // get average of ratio
 566     Float_t sum = 0;
 567     for (Int_t i=2; i<=150; ++i)
 568       sum += TMath::Abs(ratio->GetBinContent(i));
 569     sum /= 149;
 570
 571     printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
 572   }
 573
 574   TLine* line = new TLine(0, 0, displayRange, 0);
 575   line->SetLineWidth(2);
 576   line->Draw();
 577
 578   line = new TLine(0, 0.1, displayRange, 0.1);
 579   line->SetLineWidth(2);
 580   line->SetLineStyle(2);
 581   line->Draw();
 582   line = new TLine(0, -0.1, displayRange, -0.1);
 583   line->SetLineWidth(2);
 584   line->SetLineStyle(2);
 585   line->Draw();
 586
 587   canvas->Modified();
 588
 589   canvas->SaveAs(epsName);
 590   canvas->SaveAs(Form("%s.gif", epsName.Data()));
 591
 592   return canvas;
 593 }
 594
 595 void DrawResiduals(TH1* measured, TH1* unfoldedFolded, const char* epsName)
 596 {
 597   // normalize
 598   unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, 200));
 599   unfoldedFolded->Scale(measured->Integral(2, 200));
 600
 601   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
 602   canvas->Range(0, 0, 1, 1);
 603
 604   TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
 605   pad1->Draw();
 606   pad1->SetGridx();
 607   pad1->SetGridy();
 608
 609   TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
 610   pad2->Draw();
 611   pad2->SetGridx();
 612   pad2->SetGridy();
 613
 614   TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
 615   pad3->SetGridx();
 616   pad3->SetGridy();
 617   pad3->SetRightMargin(0.05);
 618   pad3->SetTopMargin(0.05);
 619   pad3->Draw();
 620
 621   pad1->SetRightMargin(0.05);
 622   pad2->SetRightMargin(0.05);
 623
 624   // no border between them
 625   pad1->SetBottomMargin(0);
 626   pad2->SetTopMargin(0);
 627
 628   pad1->cd();
 629
 630   measured->GetXaxis()->SetLabelSize(0.06);
 631   measured->GetYaxis()->SetLabelSize(0.06);
 632   measured->GetXaxis()->SetTitleSize(0.06);
 633   measured->GetYaxis()->SetTitleSize(0.06);
 634   measured->GetYaxis()->SetTitleOffset(0.6);
 635
 636   measured->GetXaxis()->SetRangeUser(0, 150);
 637
 638   measured->SetTitle(";measured multiplicity;Entries");
 639   measured->SetStats(kFALSE);
 640
 641   measured->DrawCopy("HIST");
 642   gPad->SetLogy();
 643
 644   unfoldedFolded->SetMarkerStyle(5);
 645   unfoldedFolded->SetMarkerColor(2);
 646   unfoldedFolded->SetLineColor(0);
 647   unfoldedFolded->DrawCopy("SAME P");
 648
 649   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
 650   legend->AddEntry(measured, "measured distribution");
 651   legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution");
 652   legend->SetFillColor(0);
 653   legend->Draw();
 654
 655   pad2->cd();
 656   pad2->SetBottomMargin(0.15);
 657
 658   // calculate ratio
 659   measured->Sumw2();
 660   TH1* residual = (TH1*) measured->Clone("residual");
 661   unfoldedFolded->Sumw2();
 662
 663   residual->Add(unfoldedFolded, -1);
 664
 665   // projection
 666   TH1* residualHist = new TH1F("residualHist", ";", 15, -3, 3);
 667
 668   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
 669   {
 670     if (measured->GetBinError(i) > 0)
 671     {
 672       residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
 673       residual->SetBinError(i, 1);
 674
 675       residualHist->Fill(residual->GetBinContent(i));
 676     }
 677     else
 678     {
 679       residual->SetBinContent(i, 0);
 680       residual->SetBinError(i, 0);
 681     }
 682   }
 683
 684   residual->GetYaxis()->SetTitle("Residuals   1/e (M - R #otimes U)");
 685   residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
 686   residual->DrawCopy();
 687
 688   TLine* line = new TLine(-0.5, 0, 150.5, 0);
 689   line->SetLineWidth(2);
 690   line->Draw();
 691
 692   pad3->cd();
 693   residualHist->SetStats(kFALSE);
 694   residualHist->GetXaxis()->SetLabelSize(0.08);
 695   residualHist->Fit("gaus");
 696   residualHist->Draw();
 697
 698   canvas->Modified();
 699   canvas->SaveAs(canvas->GetName());
 700
 701   //const char* epsName2 = "proj.eps";
 702   //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
 703   //canvas->SetGridx();
 704   //canvas->SetGridy();
 705
 706   //canvas->SaveAs(canvas->GetName());
 707 }
 708
 709 void bayesianExample()
 710 {
 711   TStopwatch watch;
 712   watch.Start();
 713
 714   gSystem->Load("libPWG0base");
 715
 716   TFile::Open(correctionFile);
 717   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 718   mult->LoadHistograms("Multiplicity");
 719
 720   TFile::Open(measuredFile);
 721   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 722   mult2->LoadHistograms("Multiplicity");
 723
 724   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 725
 726   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
 727
 728   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 729   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 730
 731   //mult->DrawComparison("bayesianExample", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
 732   DrawResultRatio(mcHist, result, "bayesianExample.eps");
 733
 734   //Printf("KolmogorovTest says PROB = %f", mcHist->KolmogorovTest(result, "D"));
 735   //Printf("Chi2Test says PROB = %f", mcHist->Chi2Test(result));
 736
 737   // draw residual plot
 738
 739   // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
 740   TH2* convoluted = mult->CalculateMultiplicityESD(result, etaRange);
 741   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
 742
 743   TH1* measured = mult2->GetMultiplicityESD(etaRange)->ProjectionY("measured");
 744
 745   DrawResiduals(measured, convolutedProj, "bayesianResiduals.eps");
 746
 747   watch.Stop();
 748   watch.Print();
 749 }
 750
 751 void chi2FluctuationResult()
 752 {
 753   gSystem->Load("libPWG0base");
 754
 755   TFile::Open(correctionFile);
 756   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 757   mult->LoadHistograms("Multiplicity");
 758
 759   TFile::Open(measuredFile);
 760   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 761   mult2->LoadHistograms("Multiplicity");
 762
 763   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 764   mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
 765   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 766
 767   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 768   //TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 769
 770   mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
 771
 772   TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_3");
 773   canvas->SaveAs("chi2FluctuationResult.eps");
 774 }
 775
 776 void chi2Example()
 777 {
 778   TStopwatch watch;
 779   watch.Start();
 780
 781   gSystem->Load("libPWG0base");
 782
 783   TFile::Open(correctionFile);
 784   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 785   mult->LoadHistograms("Multiplicity");
 786
 787   TFile::Open(measuredFile);
 788   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 789   mult2->LoadHistograms("Multiplicity");
 790
 791   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 792   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
 793   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 794
 795   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 796   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 797
 798   DrawResultRatio(mcHist, result, "chi2Example.eps");
 799
 800   watch.Stop();
 801   watch.Print();
 802 }
 803
 804 void chi2ExampleTPC()
 805 {
 806   TStopwatch watch;
 807   watch.Start();
 808
 809   gSystem->Load("libPWG0base");
 810
 811   TFile::Open(correctionFileTPC);
 812   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 813   mult->LoadHistograms("Multiplicity");
 814
 815   TFile::Open(measuredFileTPC);
 816   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 817   mult2->LoadHistograms("Multiplicity");
 818
 819   mult->SetMultiplicityESD(etaRangeTPC, mult2->GetMultiplicityESD(etaRangeTPC));
 820   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
 821   mult->ApplyMinuitFit(etaRangeTPC, kFALSE, AliMultiplicityCorrection::kTrVtx);
 822
 823   TH1* mcHist = mult2->GetMultiplicityVtx(etaRangeTPC)->ProjectionY("mymc");
 824   TH1* result = mult->GetMultiplicityESDCorrected(etaRangeTPC);
 825
 826   DrawResultRatio(mcHist, result, "chi2ExampleTPC.eps");
 827
 828   watch.Stop();
 829   watch.Print();
 830 }
 831
 832 void bayesianNBD()
 833 {
 834   gSystem->Load("libPWG0base");
 835   TFile::Open("multiplicityMC_3M.root");
 836
 837   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 838   mult->LoadHistograms("Multiplicity");
 839
 840   TFile::Open("multiplicityMC_3M_NBD.root");
 841   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 842   mult2->LoadHistograms("Multiplicity");
 843
 844   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 845   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1.0, 100);
 846
 847   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 848
 849   mcHist->Sumw2();
 850   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 851
 852   //mult->DrawComparison("bayesianNBD", etaRange, kFALSE, kTRUE, mcHist);
 853   DrawResultRatio(mcHist, result, "bayesianNBD.eps");
 854 }
 855
 856 void minimizationNBD()
 857 {
 858   gSystem->Load("libPWG0base");
 859   TFile::Open("multiplicityMC_3M.root");
 860
 861   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 862   mult->LoadHistograms("Multiplicity");
 863
 864   TFile::Open("multiplicityMC_3M_NBD.root");
 865   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 866   mult2->LoadHistograms("Multiplicity");
 867
 868   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 869   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
 870   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 871
 872   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 873
 874   mcHist->Sumw2();
 875   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 876
 877   //mult->DrawComparison("minimizationNBD", etaRange, kFALSE, kTRUE, mcHist);
 878   DrawResultRatio(mcHist, result, "minimizationNBD.eps");
 879 }
 880
 881 void minimizationInfluenceAlpha()
 882 {
 883   gSystem->Load("libPWG0base");
 884
 885   TFile::Open(measuredFile);
 886   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 887   mult2->LoadHistograms("Multiplicity");
 888
 889   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 890   mcHist->Scale(1.0 / mcHist->Integral());
 891   mcHist->GetXaxis()->SetRangeUser(0, 200);
 892   mcHist->SetStats(kFALSE);
 893   mcHist->SetTitle(";true multiplicity;P_{N}");
 894
 895   TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 300);
 896   canvas->Divide(3, 1);
 897
 898   TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
 899
 900   TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_100.000000");
 901   TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_03_2_100000.000000");
 902   TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_06_2_100000000.000000");
 903
 904   mcHist->Rebin(2);  mcHist->Scale(0.5);
 905   hist1->Rebin(2);   hist1->Scale(0.5);
 906   hist2->Rebin(2);   hist2->Scale(0.5);
 907   hist3->Rebin(2);   hist3->Scale(0.5);
 908
 909   mcHist->GetXaxis()->SetRangeUser(0, 200);
 910
 911   canvas->cd(1);
 912   gPad->SetLogy();
 913   mcHist->Draw();
 914   hist1->SetMarkerStyle(5);
 915   hist1->SetMarkerColor(2);
 916   hist1->Draw("SAME PE");
 917
 918   canvas->cd(2);
 919   gPad->SetLogy();
 920   mcHist->Draw();
 921   hist2->SetMarkerStyle(5);
 922   hist2->SetMarkerColor(2);
 923   hist2->Draw("SAME PE");
 924
 925   canvas->cd(3);
 926   gPad->SetLogy();
 927   mcHist->Draw();
 928   hist3->SetMarkerStyle(5);
 929   hist3->SetMarkerColor(2);
 930   hist3->Draw("SAME PE");
 931
 932   canvas->SaveAs("minimizationInfluenceAlpha.eps");
 933 }
 934
 935 void NBDFit()
 936 {
 937   gSystem->Load("libPWG0base");
 938
 939   TFile::Open(correctionFile);
 940   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 941   mult->LoadHistograms("Multiplicity");
 942
 943   TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
 944   fCurrentESD->Sumw2();
 945   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 946
 947   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])");
 948   func->SetParNames("scaling", "averagen", "k");
 949   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
 950   func->SetParLimits(1, 0.001, 1000);
 951   func->SetParLimits(2, 0.001, 1000);
 952   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
 953
 954   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
 955   lognormal->SetParNames("scaling", "mean", "sigma");
 956   lognormal->SetParameters(1, 1, 1);
 957   lognormal->SetParLimits(0, 0, 10);
 958   lognormal->SetParLimits(1, 0, 100);
 959   lognormal->SetParLimits(2, 1e-3, 10);
 960
 961   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
 962   fCurrentESD->SetStats(kFALSE);
 963   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
 964   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
 965   fCurrentESD->Draw("HIST");
 966   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
 967   fCurrentESD->Fit(func, "W0", "", 0, 50);
 968   func->SetRange(0, 100);
 969   func->Draw("SAME");
 970   printf("chi2 = %f\n", func->GetChisquare());
 971
 972   fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
 973   lognormal->SetLineColor(2);
 974   lognormal->SetLineStyle(2);
 975   lognormal->SetRange(0, 100);
 976   lognormal->Draw("SAME");
 977
 978   canvas->SaveAs("NBDFit.eps");
 979 }
 980
 981 void DifferentSamples()
 982 {
 983   // data generated by runMultiplicitySelector.C DifferentSamples
 984
 985   const char* name = "DifferentSamples";
 986
 987   TFile* file = TFile::Open(Form("%s.root", name));
 988
 989   TCanvas* canvas = new TCanvas(name, name, 800, 600);
 990   canvas->Divide(2, 2);
 991
 992   for (Int_t i=0; i<4; ++i)
 993   {
 994     canvas->cd(i+1);
 995     gPad->SetTopMargin(0.05);
 996     gPad->SetRightMargin(0.05);
 997     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", i));
 998     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", i));
 999     TH1* mc = (TH1*) file->Get(Form("mc_%d", i));
1000     mc->Sumw2();
1001
1002     chi2Result->Divide(chi2Result, mc, 1, 1, "");
1003     bayesResult->Divide(bayesResult, mc, 1, 1, "");
1004
1005     chi2Result->SetTitle(";true multiplicity;unfolded measured/MC");
1006     chi2Result->GetXaxis()->SetRangeUser(0, 150);
1007     chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1008     chi2Result->GetYaxis()->SetTitleOffset(1.2);
1009     chi2Result->SetLineColor(1);
1010     chi2Result->SetStats(kFALSE);
1011
1012     bayesResult->SetStats(kFALSE);
1013     bayesResult->SetLineColor(2);
1014
1015     chi2Result->DrawCopy("HIST");
1016     bayesResult->DrawCopy("SAME HIST");
1017
1018     TLine* line = new TLine(0, 1, 150, 1);
1019     line->Draw();
1020   }
1021
1022   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1023 }
1024
1025 void StartingConditions()
1026 {
1027   // data generated by runMultiplicitySelector.C StartingConditions
1028
1029   const char* name = "StartingConditions";
1030
1031   TFile* file = TFile::Open(Form("%s.root", name));
1032
1033   TCanvas* canvas = new TCanvas(name, name, 800, 400);
1034   canvas->Divide(2, 1);
1035
1036   TH1* mc = (TH1*) file->Get("mc");
1037   mc->Sumw2();
1038   mc->Scale(1.0 / mc->Integral());
1039
1040   //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
1041
1042   TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.95);
1043   legend->SetFillColor(0);
1044
1045   const char* names[] = { "True", "Measured 1", "Measured 2", "Measured 3", "NBD", "Flat" };
1046
1047   for (Int_t i=0; i<6; ++i)
1048   {
1049     Int_t id = i;
1050     if (id > 2)
1051       id += 2;
1052
1053     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1054     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
1055
1056     chi2Result->Divide(chi2Result, mc, 1, 1, "");
1057     bayesResult->Divide(bayesResult, mc, 1, 1, "");
1058
1059     chi2Result->SetTitle("a) #chi^{2} minimization;true multiplicity;unfolded / MC");
1060     chi2Result->GetXaxis()->SetRangeUser(0, 150);
1061     chi2Result->GetYaxis()->SetRangeUser(0.8, 1.2);
1062     chi2Result->GetYaxis()->SetTitleOffset(1.5);
1063     //chi2Result->SetMarkerStyle(marker[i]);
1064     chi2Result->SetLineColor(i+1);
1065     chi2Result->SetMarkerColor(i+1);
1066     chi2Result->SetStats(kFALSE);
1067
1068     bayesResult->SetTitle("b) Bayesian method;true multiplicity;unfolded / MC");
1069     bayesResult->GetXaxis()->SetRangeUser(0, 150);
1070     bayesResult->GetYaxis()->SetRangeUser(0.8, 1.2);
1071     bayesResult->GetYaxis()->SetTitleOffset(1.5);
1072     bayesResult->SetStats(kFALSE);
1073     //bayesResult->SetLineColor(2);
1074     bayesResult->SetLineColor(i+1);
1075
1076     canvas->cd(1);
1077     gPad->SetLeftMargin(0.12);
1078     chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1079
1080     canvas->cd(2);
1081     gPad->SetLeftMargin(0.12);
1082     bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1083
1084     //TLine* line = new TLine(0, 1, 150, 1);
1085     //line->Draw();
1086
1087     legend->AddEntry(chi2Result, names[i]);
1088   }
1089
1090   canvas->cd(1);
1091   legend->Draw();
1092   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1093 }
1094
1095 void StatisticsPlot()
1096 {
1097   const char* name = "StatisticsPlot";
1098
1099   TFile* file = TFile::Open(Form("%s.root", name));
1100
1101   TCanvas* canvas = new TCanvas(name, name, 600, 400);
1102
1103   TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
1104   fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1105   fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1106   fitResultsChi2->Draw("AP");
1107
1108   TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1109   fitResultsChi2->Fit(f);
1110
1111   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1112
1113   TH1* mc[5];
1114   TH1* result[5];
1115
1116   const char* plotname = "chi2Result";
1117
1118   name = "StatisticsPlotRatios";
1119   canvas = new TCanvas(name, name, 600, 400);
1120
1121   for (Int_t i=0; i<5; ++i)
1122   {
1123     mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1124     result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1125
1126     result[i]->SetLineColor(i+1);
1127     result[i]->Draw(((i == 0) ? "" : "SAME"));
1128   }
1129 }
1130
1131 void SystematicLowEfficiency()
1132 {
1133   gSystem->Load("libPWG0base");
1134
1135   TFile::Open(correctionFile);
1136   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1137   mult->LoadHistograms("Multiplicity");
1138
1139   // calculate result with systematic effect
1140   TFile::Open("multiplicityMC_100k_1_loweff.root");
1141   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1142   mult2->LoadHistograms("Multiplicity");
1143
1144   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1145   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1146   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1147
1148   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1149   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1150
1151   DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
1152
1153   // calculate normal result
1154   TFile::Open("multiplicityMC_100k_1.root");
1155   mult2->LoadHistograms("Multiplicity");
1156
1157   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1158   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1159
1160   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1161
1162   DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
1163
1164   Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
1165 }
1166
1167 void SystematicMisalignment()
1168 {
1169   gSystem->Load("libPWG0base");
1170
1171   TFile::Open(correctionFile);
1172   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1173   mult->LoadHistograms("Multiplicity");
1174
1175   TFile::Open("multiplicityMC_100k_fullmis.root");
1176   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1177   mult2->LoadHistograms("Multiplicity");
1178
1179   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1180   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1181   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1182
1183   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1184   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1185
1186   DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1187 }
1188
1189 void SystematicMisalignmentTPC()
1190 {
1191   gSystem->Load("libPWG0base");
1192
1193   SetTPC();
1194
1195   TFile::Open(correctionFile);
1196   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1197   mult->LoadHistograms("Multiplicity");
1198
1199   TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1200   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1201   mult2->LoadHistograms("Multiplicity");
1202
1203   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1204   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1205   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1206
1207   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1208   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1209
1210   DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1211 }
1212
1213 void EfficiencySpecies()
1214 {
1215   gSystem->Load("libPWG0base");
1216
1217   Int_t marker[] = {24, 25, 26};
1218   Int_t color[] = {1, 2, 4};
1219
1220   // SPD TPC
1221   const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
1222   Float_t etaRange[] = {0.49, 0.9};
1223   const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1224
1225   TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1226   canvas->Divide(2, 1);
1227
1228   for (Int_t loop=0; loop<2; ++loop)
1229   {
1230     Printf("%s", fileName[loop]);
1231
1232     AliCorrection* correction[4];
1233
1234     canvas->cd(loop+1);
1235
1236     gPad->SetGridx();
1237     gPad->SetGridy();
1238     gPad->SetRightMargin(0.05);
1239     //gPad->SetTopMargin(0.05);
1240
1241     TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
1242     legend->SetFillColor(0);
1243     legend->SetEntrySeparation(0.2);
1244
1245     Float_t below = 0;
1246     Float_t total = 0;
1247
1248     TFile* file = TFile::Open(fileName[loop]);
1249     if (!file)
1250     {
1251       Printf("Could not open %s", fileName[loop]);
1252       return;
1253     }
1254
1255     for (Int_t i=0; i<3; ++i)
1256     {
1257       Printf("correction %d", i);
1258
1259       TString name; name.Form("correction_%d", i);
1260       correction[i] = new AliCorrection(name, name);
1261       correction[i]->LoadHistograms();
1262
1263       TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1264       TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1265
1266       // limit vtx axis
1267       gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
1268       meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
1269
1270       // 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)
1271       /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1272         for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1273         {
1274           gene->SetBinContent(x, 0, z, 0);
1275           gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1276           meas->SetBinContent(x, 0, z, 0);
1277           meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1278         }*/
1279
1280       // limit eta axis
1281       gene->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1282       meas->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1283
1284       TH1* genePt = gene->Project3D(Form("z_%d", i));
1285       TH1* measPt = meas->Project3D(Form("z_%d", i));
1286
1287       genePt->Sumw2();
1288       measPt->Sumw2();
1289
1290       TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1291       effPt->Reset();
1292       effPt->Divide(measPt, genePt, 1, 1, "B");
1293
1294       Int_t bin = 0;
1295       for (bin=20; bin>=1; bin--)
1296       {
1297         if (effPt->GetBinContent(bin) < 0.5)
1298           break;
1299       }
1300
1301       Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1302
1303       Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1304       Printf("%.4f of the particles are below that momentum", fraction);
1305
1306       below += genePt->Integral(1, bin);
1307       total += genePt->Integral();
1308
1309       effPt->SetLineColor(color[i]);
1310       effPt->SetMarkerColor(color[i]);
1311       effPt->SetMarkerStyle(marker[i]);
1312
1313       effPt->GetXaxis()->SetRangeUser(0.06, 1);
1314       effPt->GetYaxis()->SetRangeUser(0, 1);
1315
1316       effPt->GetYaxis()->SetTitleOffset(1.2);
1317
1318       effPt->SetStats(kFALSE);
1319       effPt->SetTitle(titles[loop]);
1320       effPt->GetYaxis()->SetTitle("Efficiency");
1321
1322       effPt->DrawCopy((i == 0) ? "" : "SAME");
1323
1324       legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")));
1325     }
1326
1327     Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1328
1329     legend->Draw();
1330   }
1331
1332   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1333 }
1334
1335 void ParticleSpeciesComparison1(Bool_t chi2 = kTRUE, const char* fileNameMC = "multiplicityMC_400k_syst_species.root", const char* fileNameESD = "multiplicityMC_100k_syst.root")
1336 {
1337   gSystem->Load("libPWG0base");
1338
1339   TFile::Open(fileNameESD);
1340   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
1341   TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
1342
1343   TH1* results[10];
1344
1345   // loop over cases (normal, enhanced/reduced ratios)
1346   Int_t nMax = 7;
1347   for (Int_t i = 0; i<nMax; ++i)
1348   {
1349     TString folder;
1350     folder.Form("Multiplicity_%d", i);
1351
1352     AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
1353
1354     TFile::Open(fileNameMC);
1355     mult->LoadHistograms();
1356
1357     mult->SetMultiplicityESD(etaRange, hist);
1358
1359     if (chi2)
1360     {
1361       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1362       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1363       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1364     }
1365     else
1366     {
1367       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1368       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1369     }
1370
1371     //Float_t averageRatio = 0;
1372     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1373
1374     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1375
1376     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1377   }
1378
1379   DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
1380
1381   TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
1382
1383   for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1384   {
1385     results[0]->SetBinError(i, 0);
1386     mc->SetBinError(i, 0);
1387   }
1388
1389   const char* legendStrings[] = { "#pi^{#pm}", 0, "K^{#pm}", 0, "p,#bar{p}", 0 };
1390
1391   DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1392
1393   //not valid: draw chi2 uncertainty on top!
1394   /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1395   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1396   errorHist->SetLineColor(1);
1397   errorHist->SetLineWidth(2);
1398   TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1399   for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1400   {
1401     errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1402     errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1403   }
1404
1405   errorHist->DrawCopy("SAME");
1406   errorHist2->DrawCopy("SAME");*/
1407
1408   //canvas->SaveAs(canvas->GetName());
1409
1410   DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1411
1412   //errorHist->DrawCopy("SAME");
1413   //errorHist2->DrawCopy("SAME");
1414
1415   //canvas2->SaveAs(canvas2->GetName());
1416 }
1417
1418 /*void ParticleSpeciesComparison2()
1419 {
1420   gSystem->Load("libPWG0base");
1421
1422   const char* fileNameMC = "multiplicityMC_400k_syst.root";
1423   const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1424   Bool_t chi2 = 0;
1425
1426   TFile::Open(fileNameMC);
1427   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1428   mult->LoadHistograms();
1429
1430   TH1* mc[10];
1431   TH1* results[10];
1432
1433   // loop over cases (normal, enhanced/reduced ratios)
1434   Int_t nMax = 7;
1435   for (Int_t i = 0; i<nMax; ++i)
1436   {
1437     TString folder;
1438     folder.Form("Multiplicity_%d", i);
1439
1440     AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1441
1442     TFile::Open(fileNameESD);
1443     mult2->LoadHistograms();
1444
1445     mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1446
1447     if (chi2)
1448     {
1449       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1450       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1451       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1452     }
1453     else
1454     {
1455       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1456       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1457     }
1458
1459     //Float_t averageRatio = 0;
1460     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1461
1462     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1463
1464     TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1465     mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1466
1467     //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1468     //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1469
1470     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1471   }
1472
1473   DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1474 }*/
1475
1476 TH1* Invert(TH1* eff)
1477 {
1478   // calculate corr = 1 / eff
1479
1480   TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1481   corr->Reset();
1482
1483   for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1484   {
1485     if (eff->GetBinContent(i) > 0)
1486     {
1487       corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1488       corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1489     }
1490   }
1491
1492   return corr;
1493 }
1494
1495 void TriggerVertexCorrection()
1496 {
1497   //
1498   // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1499   //
1500
1501   gSystem->Load("libPWG0base");
1502
1503   TFile::Open(correctionFile);
1504   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1505   mult->LoadHistograms("Multiplicity");
1506
1507   TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1508   TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1509
1510   TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
1511
1512   corrINEL->SetStats(kFALSE);
1513   corrINEL->GetXaxis()->SetRangeUser(0, 30);
1514   corrINEL->GetYaxis()->SetRangeUser(0, 5);
1515   corrINEL->SetTitle(";true multiplicity;correction factor");
1516   corrINEL->SetMarkerStyle(22);
1517   corrINEL->Draw("PE");
1518
1519   corrMB->SetStats(kFALSE);
1520   corrMB->SetLineColor(2);
1521   corrMB->SetMarkerStyle(25);
1522   corrMB->SetMarkerColor(2);
1523   corrMB->Draw("SAME PE");
1524
1525   TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
1526   legend->SetFillColor(0);
1527   legend->AddEntry(corrINEL, "correction to inelastic sample");
1528   legend->AddEntry(corrMB, "correction to minimum bias sample");
1529
1530   legend->Draw();
1531
1532   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1533 }
1534
1535 void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1536 {
1537   gSystem->Load("libPWG0base");
1538
1539   TFile::Open(correctionFile);
1540   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1541   mult->LoadHistograms("Multiplicity");
1542
1543   TFile::Open(measuredFile);
1544   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1545   mult2->LoadHistograms("Multiplicity");
1546
1547   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1548
1549   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1550
1551   TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1552
1553   TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1554   TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1555
1556   if (!mc)
1557   {
1558     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1559     DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1560   }
1561
1562   TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1563   canvas->SetGridx();
1564   canvas->SetGridy();
1565   canvas->SetRightMargin(0.05);
1566   canvas->SetTopMargin(0.05);
1567
1568   errorResponse->SetLineColor(1);
1569   errorResponse->GetXaxis()->SetRangeUser(0, 200);
1570   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1571   errorResponse->SetStats(kFALSE);
1572   errorResponse->SetTitle(";true multiplicity;Uncertainty");
1573
1574   errorResponse->Draw();
1575
1576   errorMeasured->SetLineColor(2);
1577   errorMeasured->Draw("SAME");
1578
1579   errorBoth->SetLineColor(4);
1580   errorBoth->Draw("SAME");
1581
1582   Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
1583   Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
1584   Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
1585
1586   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1587
1588   TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
1589   errorResponse->Write();
1590   errorMeasured->Write();
1591   errorBoth->Write();
1592   file->Close();
1593 }
1594
1595 void StatisticalUncertaintyCompare(const char* det = "SPD")
1596 {
1597   TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1598   TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1599   TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1600   TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1601
1602   TString str;
1603   str.Form("StatisticalUncertaintyCompare%s", det);
1604
1605   TCanvas* canvas = new TCanvas(str, str, 600, 400);
1606   canvas->SetGridx();
1607   canvas->SetGridy();
1608   canvas->SetRightMargin(0.05);
1609   canvas->SetTopMargin(0.05);
1610
1611   errorResponse->SetLineColor(1);
1612   errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
1613   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1614   errorResponse->SetStats(kFALSE);
1615   errorResponse->GetYaxis()->SetTitleOffset(1.2);
1616   errorResponse->SetTitle(";true multiplicity;#sigma(U-T)/T");
1617
1618   errorResponse->Draw();
1619
1620   errorMeasured->SetLineColor(2);
1621   errorMeasured->Draw("SAME");
1622
1623   errorBoth->SetLineColor(4);
1624   errorBoth->Draw("SAME");
1625
1626   TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
1627   TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
1628
1629   errorBoth2->SetLineColor(4);
1630   errorBoth2->SetLineStyle(2);
1631   errorBoth2->Draw("SAME");
1632
1633   TLegend* legend = new TLegend(0.2, 0.6, 0.6, 0.9);
1634   legend->SetFillColor(0);
1635   legend->AddEntry(errorResponse, "response matrix (Bayesian)");
1636   legend->AddEntry(errorMeasured, "measured (Bayesian)");
1637   legend->AddEntry(errorBoth, "both (Bayesian)");
1638   legend->AddEntry(errorBoth2, "both (#chi^{2} minimization)");
1639   legend->Draw();
1640
1641   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1642 }
1643
1644 void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
1645 {
1646  const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
1647
1648   gSystem->Load("libPWG0base");
1649
1650   TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
1651   canvas->SetGridx();
1652   canvas->SetGridy();
1653   canvas->SetRightMargin(0.05);
1654   canvas->SetTopMargin(0.05);
1655
1656   AliMultiplicityCorrection* data[4];
1657   TH1* effArray[4];
1658
1659   Int_t markers[] = { 24, 25, 26, 5 };
1660   Int_t colors[] = { 1, 2, 3, 4 };
1661
1662   TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
1663   legend->SetFillColor(0);
1664
1665   TH1* effError = 0;
1666
1667   for (Int_t i=0; i<4; ++i)
1668   {
1669     TString name;
1670     name.Form("Multiplicity_%d", i);
1671
1672     TFile::Open(files[i]);
1673     data[i] = new AliMultiplicityCorrection(name, name);
1674
1675     if (i < 3)
1676     {
1677       data[i]->LoadHistograms("Multiplicity");
1678     }
1679     else
1680       data[i]->LoadHistograms("Multiplicity_0");
1681
1682     TH1* eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
1683     effArray[i] = eff;
1684
1685     eff->GetXaxis()->SetRangeUser(0, 15);
1686     eff->GetYaxis()->SetRangeUser(0, 1.1);
1687     eff->SetStats(kFALSE);
1688     eff->SetTitle(";true multiplicity;Efficiency");
1689     eff->SetLineColor(colors[i]);
1690     eff->SetMarkerColor(colors[i]);
1691     eff->SetMarkerStyle(markers[i]);
1692
1693     if (i == 3)
1694     {
1695       for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1696         eff->SetBinError(bin, 0);
1697
1698       // loop over cross section combinations
1699       for (Int_t j=1; j<7; ++j)
1700       {
1701         AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
1702         mult->LoadHistograms(Form("Multiplicity_%d", j));
1703
1704         TH1* eff2 = mult->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType);
1705
1706         for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1707         {
1708           // TODO we could also do asymmetric errors here
1709           Float_t deviation = TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
1710
1711           eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), (Double_t) deviation));
1712         }
1713       }
1714
1715       for (Int_t bin=1; bin<=20; bin++)
1716         if (eff->GetBinContent(bin) > 0)
1717           Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1718
1719       if (uncertainty) {
1720         effError = (TH1*) eff->Clone("effError");
1721         effError->Reset();
1722
1723         for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
1724           if (eff->GetBinContent(bin) > 0)
1725             effError->SetBinContent(bin, 10.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1726
1727         effError->SetLineColor(1);
1728         effError->SetMarkerStyle(1);
1729         effError->DrawCopy("SAME HIST");
1730       }
1731     }
1732
1733     eff->SetBinContent(1, 0);
1734     eff->SetBinError(1, 0);
1735
1736     canvas->cd();
1737     if (i == 0)
1738     {
1739       eff->DrawCopy("P");
1740     }
1741     else
1742       eff->DrawCopy("SAME P");
1743
1744     legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
1745   }
1746
1747   if (uncertainty)
1748     legend->AddEntry(effError, "relative syst. uncertainty #times 10");
1749
1750   legend->Draw();
1751
1752   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1753 }
1754
1755 void ModelDependencyPlot()
1756 {
1757   gSystem->Load("libPWG0base");
1758
1759   TFile::Open("multiplicityMC_3M.root");
1760   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1761   mult->LoadHistograms("Multiplicity");
1762
1763   TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
1764
1765   TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
1766   canvas->SetGridx();
1767   canvas->SetGridy();
1768   //canvas->SetRightMargin(0.05);
1769   //canvas->SetTopMargin(0.05);
1770
1771   canvas->Divide(2, 1);
1772
1773   canvas->cd(2);
1774   gPad->SetLogy();
1775
1776   Int_t selectedMult = 30;
1777   Int_t yMax = 200000;
1778
1779   TH1* full = proj->ProjectionX("full");
1780   TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
1781
1782   full->SetStats(kFALSE);
1783   full->GetXaxis()->SetRangeUser(0, 200);
1784   full->GetYaxis()->SetRangeUser(5, yMax);
1785   full->SetTitle(";multiplicity");
1786
1787   selected->SetLineColor(0);
1788   selected->SetMarkerColor(2);
1789   selected->SetMarkerStyle(7);
1790
1791   full->Draw();
1792   selected->Draw("SAME P");
1793
1794   TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
1795   legend->SetFillColor(0);
1796   legend->AddEntry(full, "true");
1797   legend->AddEntry(selected, "measured");
1798   legend->Draw();
1799
1800   TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
1801   line->SetLineWidth(2);
1802   line->Draw();
1803
1804   canvas->cd(1);
1805   gPad->SetLogy();
1806
1807   full = proj->ProjectionY("full2");
1808   selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
1809
1810   full->SetStats(kFALSE);
1811   full->GetXaxis()->SetRangeUser(0, 200);
1812   full->GetYaxis()->SetRangeUser(5, yMax);
1813   full->SetTitle(";multiplicity");
1814
1815   full->SetLineColor(0);
1816   full->SetMarkerColor(2);
1817   full->SetMarkerStyle(7);
1818
1819   full->Draw("P");
1820   selected->Draw("SAME");
1821
1822   legend->Draw();
1823
1824   line = new TLine(selectedMult, 5, selectedMult, yMax);
1825   line->SetLineWidth(2);
1826   line->Draw();
1827
1828   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1829 }
1830
1831 void SystematicpTSpectrum()
1832 {
1833   gSystem->Load("libPWG0base");
1834
1835   TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
1836   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1837   mult->LoadHistograms("Multiplicity");
1838
1839   TFile::Open("multiplicityMC_100k_syst.root");
1840   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1841   mult2->LoadHistograms("Multiplicity");
1842
1843   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1844   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1845   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1846
1847   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1848   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1849
1850   DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
1851 }
1852
1853 // to be deleted
1854 /*void covMatrix(Bool_t mc = kTRUE)
1855 {
1856   gSystem->Load("libPWG0base");
1857
1858   TFile::Open(correctionFile);
1859   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1860   mult->LoadHistograms("Multiplicity");
1861
1862   TFile::Open(measuredFile);
1863   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1864   mult2->LoadHistograms("Multiplicity");
1865
1866   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1867
1868   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1869
1870   mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
1871 }*/
1872
1873 Double_t FitPtFunc(Double_t *x, Double_t *par)
1874 {
1875   Double_t xx = x[0];
1876
1877   Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
1878   Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
1879
1880   const Float_t kTransitionWidth = 0;
1881
1882   // power law part
1883   if (xx < par[0] - kTransitionWidth)
1884   {
1885     return val1;
1886   }
1887   /*else if (xx < par[0] + kTransitionWidth)
1888   {
1889     // smooth transition
1890     Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
1891     return (1 - factor) * val1 + factor * val2;
1892   }*/
1893   else
1894   {
1895     return val2;
1896   }
1897 }
1898
1899 void FitPtNew(const char* fileName = "TruePt14TeV.root")
1900 {
1901   gSystem->Load("libANALYSIS");
1902   gSystem->Load("libPWG0base");
1903
1904   TFile::Open(fileName);
1905
1906   TH1* genePt = (TH1*) gFile->Get("fHistPt");
1907   genePt->Sumw2();
1908
1909   // normalize by bin width
1910   for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1911     genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1912
1913   genePt->GetXaxis()->SetRangeUser(0.05, 2.0);
1914
1915   genePt->Scale(1.0 / genePt->Integral());
1916
1917   TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
1918   //func->SetLineColor(2);
1919   func->SetParameters(1, -1);
1920
1921   genePt->SetMarkerStyle(25);
1922   genePt->SetTitle("");
1923   genePt->SetStats(kFALSE);
1924   genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
1925   //func->Draw("SAME");
1926
1927   genePt->Fit(func, "0", "", 0.05, 1);
1928
1929   new TCanvas;
1930   genePt->DrawCopy("P");
1931   func->SetRange(0.02, 8);
1932   func->DrawCopy("SAME");
1933   gPad->SetLogy();
1934 }
1935
1936 void FitPt(const char* fileName = "firstplots100k_truept.root")
1937 {
1938   gSystem->Load("libPWG0base");
1939
1940   TFile::Open(fileName);
1941
1942   /*
1943   // merge corrections
1944   AliCorrection* correction[4];
1945   TList list;
1946
1947   for (Int_t i=0; i<4; ++i)
1948   {
1949     Printf("correction %d", i);
1950
1951     TString name; name.Form("correction_%d", i);
1952     correction[i] = new AliCorrection(name, name);
1953     correction[i]->LoadHistograms();
1954
1955     if (i > 0)
1956       list.Add(correction[i]);
1957   }
1958
1959   correction[0]->Merge(&list);
1960
1961   TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
1962
1963   // limit vtx, eta axis
1964   gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
1965   gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
1966
1967   TH1* genePt = gene->Project3D("z");*/
1968   TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
1969   if (!genePt)
1970     genePt = (TH1*) gFile->Get("fHistPt");
1971
1972   genePt->Sumw2();
1973
1974   //genePt->Scale(1.0 / genePt->Integral());
1975
1976   // normalize by bin width
1977   for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1978     genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1979
1980   /// genePt->GetXaxis()->GetBinCenter(x));
1981
1982   genePt->GetXaxis()->SetRangeUser(0, 7.9);
1983   //genePt->GetYaxis()->SetTitle("a.u.");
1984
1985   //TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
1986   TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
1987   //func->SetLineColor(2);
1988   func->SetParameters(1, -1, 1, 1, 1);
1989   func->SetParLimits(3, 1, 10);
1990   func->SetParLimits(4, 0, 10);
1991
1992   //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
1993
1994   //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
1995   //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
1996   //func->FixParameter(0, 0.314);
1997   //func->SetParLimits(0, 0.1, 0.3);
1998
1999   genePt->SetMarkerStyle(25);
2000   genePt->SetTitle("");
2001   genePt->SetStats(kFALSE);
2002   genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
2003   //func->Draw("SAME");
2004
2005   // fit only exp. part
2006   func->SetParameters(1, -1);
2007   func->FixParameter(2, 0);
2008   func->FixParameter(3, 1);
2009   func->FixParameter(4, 1);
2010   genePt->Fit(func, "0", "", 0.2, 1);
2011
2012   new TCanvas;
2013   genePt->DrawCopy("P");
2014   func->SetRange(0.02, 8);
2015   func->DrawCopy("SAME");
2016   gPad->SetLogy();
2017
2018   // now fix exp. parameters and fit second part
2019   Double_t param0 = func->GetParameter(0);
2020   Double_t param1 = func->GetParameter(1);
2021   func->SetParameters(0, -1, 1, 1, 1);
2022   func->FixParameter(0, 0);
2023   func->FixParameter(1, -1);
2024   func->ReleaseParameter(2);
2025   func->ReleaseParameter(3);
2026   func->ReleaseParameter(4);
2027   func->SetParLimits(3, 1, 10);
2028   func->SetParLimits(4, 0, 10);
2029
2030   genePt->Fit(func, "0", "", 1.5, 4);
2031
2032   new TCanvas;
2033   genePt->DrawCopy("P");
2034   func->SetRange(0.02, 8);
2035   func->DrawCopy("SAME");
2036   gPad->SetLogy();
2037
2038   // fit both
2039   func->SetParameter(0, param0);
2040   func->SetParameter(1, param1);
2041   func->ReleaseParameter(0);
2042   func->ReleaseParameter(1);
2043
2044   new TCanvas;
2045   genePt->DrawCopy("P");
2046   func->SetRange(0.02, 5);
2047   func->DrawCopy("SAME");
2048   gPad->SetLogy();
2049
2050   genePt->Fit(func, "0", "", 0.2, 4);
2051
2052   TCanvas* canvas = new TCanvas("FitPt", "FitPt", 800, 400);
2053   canvas->Divide(2, 1);
2054   canvas->cd(1);
2055
2056   gPad->SetGridx();
2057   gPad->SetGridy();
2058   gPad->SetLeftMargin(0.13);
2059   gPad->SetRightMargin(0.05);
2060   gPad->SetTopMargin(0.05);
2061
2062   genePt->GetXaxis()->SetRangeUser(0, 4.9);
2063   genePt->GetYaxis()->SetRangeUser(1e-2, 1e4);
2064   genePt->GetYaxis()->SetTitleOffset(1.4);
2065   genePt->GetXaxis()->SetTitleOffset(1.1);
2066   genePt->DrawCopy("P");
2067   func->SetRange(0.02, 5);
2068   func->DrawCopy("SAME");
2069   gPad->SetLogy();
2070
2071   canvas->cd(2);
2072
2073   TH1* genePtClone = (TH1*) genePt->Clone("genePtClone");
2074   genePtClone->Reset();
2075   genePtClone->DrawCopy("P");
2076
2077   gPad->SetGridx();
2078   gPad->SetGridy();
2079   gPad->SetLeftMargin(0.13);
2080   gPad->SetRightMargin(0.05);
2081   gPad->SetTopMargin(0.05);
2082
2083   func->DrawCopy("SAME");
2084   gPad->SetLogy();
2085
2086   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2087
2088   TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2089
2090   TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2091
2092   TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2093
2094   for (Int_t param=0; param<5; param++)
2095   {
2096     for (Int_t sign=0; sign<2; sign++)
2097     {
2098       TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign));  // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2099       func2->SetParameters(func->GetParameters());
2100       //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2101
2102       Float_t factor = ((sign == 0) ? 0.9 : 1.1);
2103       func2->SetParameter(param, func2->GetParameter(param) * factor);
2104       //func2->Print();
2105
2106       canvas->cd(2);
2107       func2->SetLineWidth(1);
2108       func2->SetLineColor(2);
2109       func2->DrawCopy("SAME");
2110
2111       canvas2->cd();
2112       TH1* second = func2->GetHistogram();
2113       second->Divide(first);
2114       second->SetLineColor(param + 1);
2115       second->GetYaxis()->SetRangeUser(0, 2);
2116       second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2117       second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2118     }
2119   }
2120
2121   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2122   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2123
2124   file->Close();
2125 }
2126
2127 void DrawSystematicpT()
2128 {
2129   TFile* file = TFile::Open("SystematicpT.root");
2130
2131   TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2132   TH1* result2 = (TH1*) file->Get("result_unity");
2133
2134   TH1* mcHist[12];
2135   TH1* result[12];
2136
2137   Int_t nParams = 5;
2138
2139   for (Int_t id=0; id<nParams*2; ++id)
2140   {
2141     mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2142     result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2143   }
2144
2145   DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2146
2147   //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2148
2149   DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2150
2151   //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2152
2153   // does not make sense: mc is different
2154   //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2155 }
2156
2157 void SystematicpT(Bool_t chi2 = 1)
2158 {
2159   gSystem->Load("libPWG0base");
2160
2161   TFile::Open("ptspectrum900.root");
2162   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2163   mult->LoadHistograms("Multiplicity");
2164
2165   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2166
2167   TH1* mcHist[12];
2168   TH1* result[12];
2169
2170   Int_t nParams = 5;
2171
2172   for (Int_t param=0; param<nParams; param++)
2173   {
2174     for (Int_t sign=0; sign<2; sign++)
2175     {
2176       // calculate result with systematic effect
2177       TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2178       mult2->LoadHistograms("Multiplicity");
2179
2180       mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2181
2182       if (chi2)
2183       {
2184         mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2185         mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2186       }
2187       else
2188         mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2189
2190       Int_t id = param * 2 + sign;
2191
2192       mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2193       result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2194
2195       TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2196       DrawResultRatio(mcHist[id], result[id], tmp);
2197     }
2198   }
2199
2200   // calculate normal result
2201   TFile::Open("ptspectrum100_1.root");
2202   mult2->LoadHistograms("Multiplicity");
2203   TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2204
2205   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2206
2207   if (chi2)
2208   {
2209     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2210   }
2211   else
2212     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2213
2214   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2215
2216   TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2217   mcHist2->Write();
2218   result2->Write();
2219   for (Int_t id=0; id<nParams*2; ++id)
2220   {
2221     mcHist[id]->Write();
2222     result[id]->Write();
2223   }
2224   file->Close();
2225
2226   DrawSystematicpT();
2227 }
2228
2229 void DrawSystematicpT2()
2230 {
2231   //displayRange = 200;
2232
2233   // read from file
2234   TFile* file = TFile::Open("SystematicpT2.root");
2235   TH1* mcHist = (TH1*) file->Get("mymc");
2236   TH1* result[12];
2237   result[0] = (TH1*) file->Get("result_unity");
2238   Int_t nParams = 5;
2239   for (Int_t id=0; id<nParams*2; ++id)
2240     result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2241
2242   DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2243   DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2244   DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2245 }
2246
2247 void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2248 {
2249   gSystem->Load("libPWG0base");
2250
2251   if (tpc)
2252   {
2253     SetTPC();
2254     TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2255   }
2256   else
2257     TFile::Open("ptspectrum100_1.root");
2258
2259   AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2260   measured->LoadHistograms("Multiplicity");
2261   TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2262
2263   TH1* result[12];
2264
2265   Int_t nParams = 5;
2266
2267   // -1 = unity change, 0...4 parameters
2268   for (Int_t id=-1; id<nParams*2; id++)
2269   {
2270     Int_t param = id / 2;
2271     Int_t sign = id % 2;
2272
2273     TString idStr;
2274     if (id == -1)
2275     {
2276       idStr = "unity";
2277     }
2278     else
2279       idStr.Form("%d_%d", param, sign);
2280
2281     // calculate result with systematic effect
2282     if (tpc)
2283     {
2284       TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2285     }
2286     else
2287       TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2288
2289     AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2290     response->LoadHistograms("Multiplicity");
2291
2292     response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2293
2294     if (chi2)
2295     {
2296       response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2297       response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2298     }
2299     else
2300       response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2301
2302     result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2303
2304     TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2305     DrawResultRatio(mcHist, result[id+1], tmp);
2306   }
2307
2308   TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2309   mcHist->Write();
2310   for (Int_t id=0; id<nParams*2+1; ++id)
2311     result[id]->Write();
2312   file->Close();
2313
2314   DrawSystematicpT2();
2315 }
2316
2317 void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2318 {
2319   // only needed for TPC
2320   SetTPC();
2321
2322   gSystem->Load("libPWG0base");
2323
2324   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2325   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2326   mult->LoadHistograms("Multiplicity");
2327
2328   TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2329   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2330   mult2->LoadHistograms("Multiplicity");
2331
2332   // "normal" result
2333   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2334
2335   if (chi2)
2336   {
2337     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2338     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2339   }
2340   else
2341     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2342
2343   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2344   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2345
2346   TH1* syst[2];
2347
2348   // change of pt spectrum (down)
2349   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2350   AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2351   mult3->LoadHistograms("Multiplicity");
2352   mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2353   if (chi2)
2354   {
2355     mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2356     mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2357   }
2358   else
2359     mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2360   syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2361
2362   // change of pt spectrum (up)
2363   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2364   AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2365   mult4->LoadHistograms("Multiplicity");
2366   mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2367   if (chi2)
2368   {
2369     mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2370     mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2371   }
2372   else
2373     mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2374   syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
2375
2376   DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
2377
2378   Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2379   Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2380 }
2381
2382 TH1* SystematicsSummary(Bool_t tpc = 1)
2383 {
2384   Int_t nEffects = 7;
2385
2386   TH1* effects[10];
2387   const char** names = 0;
2388   Int_t colors[] = { 1, 2, 3, 4, 6, 7, 8 };
2389   Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
2390
2391   for (Int_t i=0; i<nEffects; ++i)
2392     effects[i] = new TH1F("SystematicsSummary", ";true multiplicity;Effect", 201, -0.5, 200.5);
2393
2394   if (tpc)
2395   {
2396     SetTPC();
2397
2398     const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2399     names = namesTPC;
2400
2401     // method
2402     TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2403     TH1* hist = (TH1*) file->Get("errorBoth");
2404
2405     // smooth a bit, but skip 0 bin
2406     effects[0]->SetBinContent(2, hist->GetBinContent(2));
2407     for (Int_t i=3; i<=200; ++i)
2408       effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2409
2410     // relative x-section
2411     effects[1]->SetBinContent(2, 0.005);
2412     effects[1]->SetBinContent(3, 0.0025);
2413     effects[1]->SetBinContent(4, 0.0025);
2414
2415     // particle composition
2416     for (Int_t i=2; i<=101; ++i)
2417     {
2418       if (i < 41)
2419       {
2420         effects[2]->SetBinContent(i, 0.01);
2421       }
2422       else if (i < 76)
2423       {
2424         effects[2]->SetBinContent(i, 0.02);
2425       }
2426       else
2427         effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2428     }
2429
2430     // pt cut off (only tpc)
2431     for (Int_t i=2; i<=101; ++i)
2432     {
2433       if (i < 11)
2434       {
2435         effects[3]->SetBinContent(i, 0.05);
2436       }
2437       else if (i < 51)
2438       {
2439         effects[3]->SetBinContent(i, 0.01);
2440       }
2441       else
2442         effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2443     }
2444
2445     // track selection (only tpc)
2446     for (Int_t i=2; i<=101; ++i)
2447       effects[4]->SetBinContent(i, 0.03);
2448
2449     // secondaries
2450     for (Int_t i=2; i<=101; ++i)
2451       effects[5]->SetBinContent(i, 0.01);
2452
2453     // pt spectrum
2454     for (Int_t i=2; i<=101; ++i)
2455     {
2456       if (i < 21)
2457       {
2458         effects[6]->SetBinContent(i, 0.05);
2459       }
2460       else if (i < 51)
2461       {
2462         effects[6]->SetBinContent(i, 0.02);
2463       }
2464       else
2465         effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2466     }
2467
2468   }
2469   else
2470   {
2471     displayRange = 200;
2472     nEffects = 5;
2473
2474     const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Rel. cross-section", "Particle composition", "Secondaries", "p_{t} spectrum"};
2475     names = namesSPD;
2476
2477     // method
2478     TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2479     TH1* hist = (TH1*) file->Get("errorBoth");
2480
2481     // smooth a bit, but skip 0 bin
2482     effects[0]->SetBinContent(2, hist->GetBinContent(2));
2483     for (Int_t i=3; i<=201; ++i)
2484       effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2485
2486     // relative x-section
2487     effects[1]->SetBinContent(2, 0.01);
2488     effects[1]->SetBinContent(3, 0.005);
2489
2490     // particle composition
2491     for (Int_t i=2; i<=201; ++i)
2492     {
2493       if (i < 6)
2494       {
2495         effects[2]->SetBinContent(i, 0.3);
2496       }
2497       else if (i < 11)
2498       {
2499         effects[2]->SetBinContent(i, 0.05);
2500       }
2501       else if (i < 121)
2502       {
2503         effects[2]->SetBinContent(i, 0.02);
2504       }
2505       else if (i < 151)
2506       {
2507         effects[2]->SetBinContent(i, 0.02 + 0.04 / 30 * (i - 121));
2508       }
2509       else
2510         effects[2]->SetBinContent(i, 0.06 + 0.1 / 30 * (i - 151));
2511     }
2512
2513     // secondaries
2514     for (Int_t i=2; i<=201; ++i)
2515       effects[3]->SetBinContent(i, 0.01);
2516
2517     // pt spectrum
2518     for (Int_t i=2; i<=201; ++i)
2519     {
2520       if (i < 6)
2521       {
2522         effects[4]->SetBinContent(i, 1);
2523       }
2524       else if (i < 121)
2525       {
2526         effects[4]->SetBinContent(i, 0.03);
2527       }
2528       else if (i < 151)
2529       {
2530         effects[4]->SetBinContent(i, 0.03 + 0.07 / 30 * (i - 121));
2531       }
2532       else
2533         effects[4]->SetBinContent(i, 0.1);
2534     }
2535   }
2536
2537   TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 400);
2538   canvas->SetRightMargin(0.25);
2539   canvas->SetTopMargin(0.05);
2540   TLegend* legend = new TLegend(0.2, 0.4, 0.5, 0.4 + 0.5 * nEffects / 7);
2541   legend->SetFillColor(0);
2542
2543   for (Int_t i=0; i<nEffects; ++i)
2544   {
2545     TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2546     /*current->Reset();
2547     for (Int_t j=0; j<nEffects-i; ++j)
2548       current->Add(effects[j]);*/
2549
2550     current->SetLineColor(colors[i]);
2551     //current->SetFillColor(colors[i]);
2552     current->SetMarkerColor(colors[i]);
2553     //current->SetMarkerStyle(markers[i]);
2554
2555     current->SetStats(kFALSE);
2556     current->GetYaxis()->SetRangeUser(0, 0.4);
2557     current->GetXaxis()->SetRangeUser(0, displayRange);
2558     current->DrawCopy(((i == 0) ? "" : "SAME"));
2559     legend->AddEntry(current, names[i]);
2560
2561     TLatex* text = new TLatex(displayRange+5, current->GetBinContent(displayRange+1), names[i]);
2562     text->SetTextColor(colors[i]);
2563     text->Draw();
2564   }
2565
2566   // add total in square
2567   TH1* total = (TH1*) effects[0]->Clone("total");
2568   total->Reset();
2569
2570   for (Int_t i=0; i<nEffects; ++i)
2571   {
2572     //Printf("%d %f", i, effects[i]->GetBinContent(20));
2573     effects[i]->Multiply(effects[i]);
2574     total->Add(effects[i]);
2575   }
2576
2577   for (Int_t i=1; i<=total->GetNbinsX(); ++i)
2578     if (total->GetBinContent(i) > 0)
2579       total->SetBinContent(i, TMath::Min(sqrt(total->GetBinContent(i)), 1.0));
2580
2581   //Printf("%f", total->GetBinContent(20));
2582
2583   total->SetMarkerStyle(3);
2584   total->SetMarkerColor(1);
2585   legend->AddEntry(total, "total");
2586   total->DrawCopy("SAME P");
2587
2588   legend->Draw();
2589
2590   canvas->SaveAs(canvas->GetName());
2591
2592   return total;
2593 }
2594
2595 void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE, Bool_t small = kFALSE)
2596 {
2597   gSystem->Load("libPWG0base");
2598
2599   if (tpc)
2600     SetTPC();
2601
2602   if (!chi2)
2603     Printf("WARNING: Bayesian set. This is only for test!");
2604
2605   // systematic error
2606   TH1* error = SystematicsSummary(tpc);
2607
2608   TFile::Open(correctionFile);
2609   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2610   mult->LoadHistograms("Multiplicity");
2611
2612   TFile::Open(measuredFile);
2613   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2614   mult2->LoadHistograms("Multiplicity");
2615
2616   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2617
2618   if (chi2)
2619   {
2620     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2621     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kINEL);
2622   }
2623   else
2624     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kINEL, 1, 100, 0, kFALSE);
2625
2626   TH1* mcHist = mult2->GetMultiplicityINEL(etaRange)->ProjectionY("mymc");
2627   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2628
2629   DrawResultRatio(mcHist, result, "finalPlotCheck.eps");
2630
2631   // normalize result
2632   result->Scale(1.0 / result->Integral(2, 200));
2633
2634   result->GetXaxis()->SetRangeUser(0, ((tpc) ? displayRange : 200));
2635   result->SetBinContent(1, 0); result->SetBinError(1, 0);
2636   result->SetTitle(";true multiplicity;Probability");
2637   result->SetLineColor(1);
2638   result->SetStats(kFALSE);
2639
2640   TH1* systError = (TH1*) result->Clone("systError");
2641   for (Int_t i=2; i<=systError->GetNbinsX(); ++i)
2642     systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2643
2644   // change error drawing style
2645   systError->SetFillColor(15);
2646
2647   TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", (small) ? 600 : 800, 400);
2648   canvas->SetRightMargin(0.05);
2649   canvas->SetTopMargin(0.05);
2650
2651   systError->Draw("E2 ][");
2652   result->DrawCopy("SAME E ][");
2653   canvas->SetLogy();
2654
2655   //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
2656   TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
2657   text->SetFillColor(0);
2658   text->SetTextAlign(12);
2659   text->AddText("Systematic errors summed quadratically");
2660   text->AddText("0.6 million minimum bias events");
2661   text->AddText("corrected to inelastic events");
2662   text->Draw("B");
2663
2664   TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
2665   text2->SetFillColor(0);
2666   text2->SetTextAlign(12);
2667   text2->AddText("#sqrt{s} = 14 TeV");
2668   if (tpc)
2669   {
2670     text2->AddText("|#eta| < 0.9");
2671   }
2672   else
2673     text2->AddText("|#eta| < 2.0");
2674   text2->AddText("simulated data (PYTHIA)");
2675   text2->Draw("B");
2676
2677   if (tpc)
2678   {
2679     TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
2680     text3->SetNDC();
2681     text3->Draw();
2682   }
2683   else
2684   {
2685     TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
2686     text3->SetNDC();
2687     text3->Draw();
2688   }
2689
2690   // alice logo
2691   TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
2692   pad->Draw();
2693   pad->cd();
2694   TImage* img = TImage::Open("$HOME/alice.png");
2695   img->SetImageQuality(TAttImage::kImgBest);
2696   img->Draw();
2697
2698   canvas->Modified();
2699
2700 /*  TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
2701   text->SetTextSize(0.04);
2702   text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
2703   text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
2704   text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
2705   text->Draw();*/
2706
2707
2708   canvas->SaveAs(canvas->GetName());
2709 }
2710
2711 void BlobelUnfoldingExample()
2712 {
2713   const Int_t kSize = 20;
2714
2715   TMatrixD matrix(kSize, kSize);
2716   for (Int_t x=0; x<kSize; x++)
2717   {
2718     for (Int_t y=0; y<kSize; y++)
2719     {
2720       if (x == y)
2721       {
2722         if (x == 0 || x == kSize -1)
2723         {
2724           matrix(x, y) = 0.75;
2725         }
2726         else
2727           matrix(x, y) = 0.5;
2728       }
2729       else if (TMath::Abs(x - y) == 1)
2730       {
2731         matrix(x, y) = 0.25;
2732       }
2733     }
2734   }
2735
2736   //matrix.Print();
2737
2738   TMatrixD inverted(matrix);
2739   inverted.Invert();
2740
2741   //inverted.Print();
2742
2743   TH1F* inputDist = new TH1F("inputDist", ";t;#tilde{T}(t)", kSize, -0.5, (Float_t) kSize - 0.5);
2744   TVectorD inputDistVector(kSize);
2745   TH1F* unfolded = inputDist->Clone("unfolded");
2746   TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
2747   measuredIdealDist->SetTitle(";m;#tilde{M}(m)");
2748   TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
2749
2750   TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
2751   // norm: 1/(sqrt(2pi)sigma)
2752   gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
2753   //gaus->Print();
2754
2755   for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
2756   {
2757     Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
2758     inputDist->SetBinContent(x, value);
2759     inputDistVector(x-1) = value;
2760   }
2761
2762   TVectorD measuredDistIdealVector = matrix * inputDistVector;
2763
2764   for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
2765     measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
2766
2767   measuredDist->FillRandom(measuredIdealDist, 10000);
2768
2769   // fill error matrix before scaling
2770   TMatrixD covarianceMatrixMeasured(kSize, kSize);
2771   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
2772     covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
2773
2774   TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
2775   //covarianceMatrix.Print();
2776
2777   TVectorD measuredDistVector(kSize);
2778   for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
2779     measuredDistVector(x-1) = measuredDist->GetBinContent(x);
2780
2781   TVectorD unfoldedVector = inverted * measuredDistVector;
2782   for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
2783     unfolded->SetBinContent(x, unfoldedVector(x-1));
2784
2785   TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1000, 500);
2786   canvas->SetTopMargin(0.05);
2787   canvas->Divide(2, 1);
2788
2789   canvas->cd(1);
2790   canvas->cd(1)->SetLeftMargin(0.15);
2791   canvas->cd(1)->SetRightMargin(0.05);
2792   measuredDist->GetYaxis()->SetTitleOffset(1.7);
2793   measuredDist->SetStats(0);
2794   measuredDist->DrawCopy();
2795   gaus->Draw("SAME");
2796
2797   canvas->cd(2);
2798   canvas->cd(2)->SetLeftMargin(0.15);
2799   canvas->cd(2)->SetRightMargin(0.05);
2800   unfolded->GetYaxis()->SetTitleOffset(1.7);
2801   unfolded->SetStats(0);
2802   unfolded->DrawCopy();
2803   gaus->Draw("SAME");
2804
2805   canvas->SaveAs("BlobelUnfoldingExample.eps");
2806 }
2807
2808 void E735Fit()
2809 {
2810   TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
2811   fCurrentESD->Sumw2();
2812
2813   // Open the input stream
2814   ifstream in;
2815   in.open("e735data.txt");
2816
2817   while(in.good())
2818   {
2819     Float_t x, y, ye;
2820     in >> x >> y >> ye;
2821
2822     //Printf("%f %f %f", x, y, ye);
2823     fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
2824     fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
2825   }
2826
2827   in.close();
2828
2829   //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
2830
2831   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2832
2833   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])");
2834   func->SetParNames("scaling", "averagen", "k");
2835   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
2836   func->SetParLimits(1, 0.001, 1000);
2837   func->SetParLimits(2, 0.001, 1000);
2838   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
2839
2840   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
2841   lognormal->SetParNames("scaling", "mean", "sigma");
2842   lognormal->SetParameters(1, 1, 1);
2843   lognormal->SetParLimits(0, 0, 10);
2844   lognormal->SetParLimits(1, 0, 100);
2845   lognormal->SetParLimits(2, 1e-3, 10);
2846
2847   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
2848   fCurrentESD->SetStats(kFALSE);
2849   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
2850   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
2851   fCurrentESD->Draw("");
2852   fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
2853   fCurrentESD->Fit(func, "0", "", 0, 150);
2854   func->SetRange(0, 250);
2855   func->Draw("SAME");
2856   printf("chi2 = %f\n", func->GetChisquare());
2857
2858   fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
2859   lognormal->SetLineColor(2);
2860   lognormal->SetLineStyle(2);
2861   lognormal->SetRange(0, 250);
2862   lognormal->Draw("SAME");
2863
2864   gPad->SetLogy();
2865
2866   canvas->SaveAs("E735Fit.eps");
2867 }