PWG0/dNdEta/plotsMultiplicity.C

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