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 = 150; // 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   // draw residual plot
 733
 734   // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
 735   TH2* convoluted = mult->CalculateMultiplicityESD(result, etaRange);
 736   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
 737
 738   TH1* measured = mult2->GetMultiplicityESD(etaRange)->ProjectionY("measured");
 739
 740   DrawResiduals(measured, convolutedProj, "bayesianResiduals.eps");
 741
 742   watch.Stop();
 743   watch.Print();
 744 }
 745
 746 void chi2FluctuationResult()
 747 {
 748   gSystem->Load("libPWG0base");
 749
 750   TFile::Open(correctionFile);
 751   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 752   mult->LoadHistograms("Multiplicity");
 753
 754   TFile::Open(measuredFile);
 755   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 756   mult2->LoadHistograms("Multiplicity");
 757
 758   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 759   mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
 760   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 761
 762   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 763   //TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 764
 765   mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
 766
 767   TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_3");
 768   canvas->SaveAs("chi2FluctuationResult.eps");
 769 }
 770
 771 void chi2Example()
 772 {
 773   TStopwatch watch;
 774   watch.Start();
 775
 776   gSystem->Load("libPWG0base");
 777
 778   TFile::Open(correctionFile);
 779   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 780   mult->LoadHistograms("Multiplicity");
 781
 782   TFile::Open(measuredFile);
 783   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 784   mult2->LoadHistograms("Multiplicity");
 785
 786   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 787   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
 788   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 789
 790   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 791   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 792
 793   DrawResultRatio(mcHist, result, "chi2Example.eps");
 794
 795   watch.Stop();
 796   watch.Print();
 797 }
 798
 799 void chi2ExampleTPC()
 800 {
 801   TStopwatch watch;
 802   watch.Start();
 803
 804   gSystem->Load("libPWG0base");
 805
 806   TFile::Open(correctionFileTPC);
 807   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 808   mult->LoadHistograms("Multiplicity");
 809
 810   TFile::Open(measuredFileTPC);
 811   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 812   mult2->LoadHistograms("Multiplicity");
 813
 814   mult->SetMultiplicityESD(etaRangeTPC, mult2->GetMultiplicityESD(etaRangeTPC));
 815   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
 816   mult->ApplyMinuitFit(etaRangeTPC, kFALSE, AliMultiplicityCorrection::kTrVtx);
 817
 818   TH1* mcHist = mult2->GetMultiplicityVtx(etaRangeTPC)->ProjectionY("mymc");
 819   TH1* result = mult->GetMultiplicityESDCorrected(etaRangeTPC);
 820
 821   DrawResultRatio(mcHist, result, "chi2ExampleTPC.eps");
 822
 823   watch.Stop();
 824   watch.Print();
 825 }
 826
 827 void bayesianNBD()
 828 {
 829   gSystem->Load("libPWG0base");
 830   TFile::Open("multiplicityMC_3M.root");
 831
 832   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 833   mult->LoadHistograms("Multiplicity");
 834
 835   TFile::Open("multiplicityMC_3M_NBD.root");
 836   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 837   mult2->LoadHistograms("Multiplicity");
 838
 839   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 840   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1.0, 100);
 841
 842   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 843
 844   mcHist->Sumw2();
 845   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 846
 847   //mult->DrawComparison("bayesianNBD", etaRange, kFALSE, kTRUE, mcHist);
 848   DrawResultRatio(mcHist, result, "bayesianNBD.eps");
 849 }
 850
 851 void minimizationNBD()
 852 {
 853   gSystem->Load("libPWG0base");
 854   TFile::Open("multiplicityMC_3M.root");
 855
 856   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 857   mult->LoadHistograms("Multiplicity");
 858
 859   TFile::Open("multiplicityMC_3M_NBD.root");
 860   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 861   mult2->LoadHistograms("Multiplicity");
 862
 863   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 864   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
 865   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 866
 867   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 868
 869   mcHist->Sumw2();
 870   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 871
 872   //mult->DrawComparison("minimizationNBD", etaRange, kFALSE, kTRUE, mcHist);
 873   DrawResultRatio(mcHist, result, "minimizationNBD.eps");
 874 }
 875
 876 void minimizationInfluenceAlpha()
 877 {
 878   gSystem->Load("libPWG0base");
 879
 880   TFile::Open(measuredFile);
 881   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 882   mult2->LoadHistograms("Multiplicity");
 883
 884   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 885   mcHist->Scale(1.0 / mcHist->Integral());
 886   mcHist->GetXaxis()->SetRangeUser(0, 200);
 887   mcHist->SetStats(kFALSE);
 888   mcHist->SetTitle(";true multiplicity;P_{N}");
 889
 890   TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 300);
 891   canvas->Divide(3, 1);
 892
 893   TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
 894
 895   TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_100.000000");
 896   TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_03_2_100000.000000");
 897   TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_06_2_100000000.000000");
 898
 899   mcHist->Rebin(2);  mcHist->Scale(0.5);
 900   hist1->Rebin(2);   hist1->Scale(0.5);
 901   hist2->Rebin(2);   hist2->Scale(0.5);
 902   hist3->Rebin(2);   hist3->Scale(0.5);
 903
 904   mcHist->GetXaxis()->SetRangeUser(0, 200);
 905
 906   canvas->cd(1);
 907   gPad->SetLogy();
 908   mcHist->Draw();
 909   hist1->SetMarkerStyle(5);
 910   hist1->SetMarkerColor(2);
 911   hist1->Draw("SAME PE");
 912
 913   canvas->cd(2);
 914   gPad->SetLogy();
 915   mcHist->Draw();
 916   hist2->SetMarkerStyle(5);
 917   hist2->SetMarkerColor(2);
 918   hist2->Draw("SAME PE");
 919
 920   canvas->cd(3);
 921   gPad->SetLogy();
 922   mcHist->Draw();
 923   hist3->SetMarkerStyle(5);
 924   hist3->SetMarkerColor(2);
 925   hist3->Draw("SAME PE");
 926
 927   canvas->SaveAs("minimizationInfluenceAlpha.eps");
 928 }
 929
 930 void NBDFit()
 931 {
 932   gSystem->Load("libPWG0base");
 933
 934   TFile::Open(correctionFile);
 935   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 936   mult->LoadHistograms("Multiplicity");
 937
 938   TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
 939   fCurrentESD->Sumw2();
 940   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 941
 942   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])");
 943   func->SetParNames("scaling", "averagen", "k");
 944   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
 945   func->SetParLimits(1, 0.001, 1000);
 946   func->SetParLimits(2, 0.001, 1000);
 947   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
 948
 949   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
 950   lognormal->SetParNames("scaling", "mean", "sigma");
 951   lognormal->SetParameters(1, 1, 1);
 952   lognormal->SetParLimits(0, 0, 10);
 953   lognormal->SetParLimits(1, 0, 100);
 954   lognormal->SetParLimits(2, 1e-3, 10);
 955
 956   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
 957   fCurrentESD->SetStats(kFALSE);
 958   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
 959   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
 960   fCurrentESD->Draw("HIST");
 961   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
 962   fCurrentESD->Fit(func, "W0", "", 0, 50);
 963   func->SetRange(0, 100);
 964   func->Draw("SAME");
 965   printf("chi2 = %f\n", func->GetChisquare());
 966
 967   fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
 968   lognormal->SetLineColor(2);
 969   lognormal->SetLineStyle(2);
 970   lognormal->SetRange(0, 100);
 971   lognormal->Draw("SAME");
 972
 973   canvas->SaveAs("NBDFit.eps");
 974 }
 975
 976 void DifferentSamples()
 977 {
 978   // data generated by runMultiplicitySelector.C DifferentSamples
 979
 980   const char* name = "DifferentSamples";
 981
 982   TFile* file = TFile::Open(Form("%s.root", name));
 983
 984   TCanvas* canvas = new TCanvas(name, name, 800, 600);
 985   canvas->Divide(2, 2);
 986
 987   for (Int_t i=0; i<4; ++i)
 988   {
 989     canvas->cd(i+1);
 990     gPad->SetTopMargin(0.05);
 991     gPad->SetRightMargin(0.05);
 992     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", i));
 993     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", i));
 994     TH1* mc = (TH1*) file->Get(Form("mc_%d", i));
 995     mc->Sumw2();
 996
 997     chi2Result->Divide(chi2Result, mc, 1, 1, "");
 998     bayesResult->Divide(bayesResult, mc, 1, 1, "");
 999
1000     chi2Result->SetTitle(";true multiplicity;unfolded measured/MC");
1001     chi2Result->GetXaxis()->SetRangeUser(0, 150);
1002     chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1003     chi2Result->GetYaxis()->SetTitleOffset(1.2);
1004     chi2Result->SetLineColor(1);
1005     chi2Result->SetStats(kFALSE);
1006
1007     bayesResult->SetStats(kFALSE);
1008     bayesResult->SetLineColor(2);
1009
1010     chi2Result->DrawCopy("HIST");
1011     bayesResult->DrawCopy("SAME HIST");
1012
1013     TLine* line = new TLine(0, 1, 150, 1);
1014     line->Draw();
1015   }
1016
1017   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1018 }
1019
1020 void StartingConditions()
1021 {
1022   // data generated by runMultiplicitySelector.C StartingConditions
1023
1024   const char* name = "StartingConditions";
1025
1026   TFile* file = TFile::Open(Form("%s.root", name));
1027
1028   TCanvas* canvas = new TCanvas(name, name, 800, 400);
1029   canvas->Divide(2, 1);
1030
1031   TH1* mc = (TH1*) file->Get("mc");
1032   mc->Sumw2();
1033   mc->Scale(1.0 / mc->Integral());
1034
1035   //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
1036
1037   for (Int_t i=0; i<6; ++i)
1038   {
1039     Int_t id = i;
1040     if (id > 2)
1041       id += 2;
1042
1043     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
1044     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
1045
1046     chi2Result->Divide(chi2Result, mc, 1, 1, "");
1047     bayesResult->Divide(bayesResult, mc, 1, 1, "");
1048
1049     chi2Result->SetTitle("a) #chi^{2} minimization;true multiplicity;unfolded / MC");
1050     chi2Result->GetXaxis()->SetRangeUser(0, 150);
1051     chi2Result->GetYaxis()->SetRangeUser(0.8, 1.2);
1052     chi2Result->GetYaxis()->SetTitleOffset(1.5);
1053     //chi2Result->SetMarkerStyle(marker[i]);
1054     chi2Result->SetLineColor(i+1);
1055     chi2Result->SetStats(kFALSE);
1056
1057     bayesResult->SetTitle("b) Bayesian method;true multiplicity;unfolded / MC");
1058     bayesResult->GetXaxis()->SetRangeUser(0, 150);
1059     bayesResult->GetYaxis()->SetRangeUser(0.8, 1.2);
1060     bayesResult->GetYaxis()->SetTitleOffset(1.5);
1061     bayesResult->SetStats(kFALSE);
1062     bayesResult->SetLineColor(2);
1063     bayesResult->SetLineColor(i+1);
1064
1065     canvas->cd(1);
1066     gPad->SetLeftMargin(0.12);
1067     chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1068
1069     canvas->cd(2);
1070     gPad->SetLeftMargin(0.12);
1071     bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
1072
1073     //TLine* line = new TLine(0, 1, 150, 1);
1074     //line->Draw();
1075   }
1076
1077   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1078 }
1079
1080 void StatisticsPlot()
1081 {
1082   const char* name = "StatisticsPlot";
1083
1084   TFile* file = TFile::Open(Form("%s.root", name));
1085
1086   TCanvas* canvas = new TCanvas(name, name, 600, 400);
1087
1088   TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
1089   fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1090   fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1091   fitResultsChi2->Draw("AP");
1092
1093   TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1094   fitResultsChi2->Fit(f);
1095
1096   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1097
1098   TH1* mc[5];
1099   TH1* result[5];
1100
1101   const char* plotname = "chi2Result";
1102
1103   name = "StatisticsPlotRatios";
1104   canvas = new TCanvas(name, name, 600, 400);
1105
1106   for (Int_t i=0; i<5; ++i)
1107   {
1108     mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1109     result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1110
1111     result[i]->SetLineColor(i+1);
1112     result[i]->Draw(((i == 0) ? "" : "SAME"));
1113   }
1114 }
1115
1116 void SystematicLowEfficiency()
1117 {
1118   gSystem->Load("libPWG0base");
1119
1120   TFile::Open(correctionFile);
1121   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1122   mult->LoadHistograms("Multiplicity");
1123
1124   // calculate result with systematic effect
1125   TFile::Open("multiplicityMC_100k_1_loweff.root");
1126   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1127   mult2->LoadHistograms("Multiplicity");
1128
1129   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1130   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1131   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1132
1133   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1134   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1135
1136   DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
1137
1138   // calculate normal result
1139   TFile::Open("multiplicityMC_100k_1.root");
1140   mult2->LoadHistograms("Multiplicity");
1141
1142   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1143   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1144
1145   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1146
1147   DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
1148
1149   Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
1150 }
1151
1152 void SystematicMisalignment()
1153 {
1154   gSystem->Load("libPWG0base");
1155
1156   TFile::Open(correctionFile);
1157   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1158   mult->LoadHistograms("Multiplicity");
1159
1160   TFile::Open("multiplicityMC_100k_fullmis.root");
1161   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1162   mult2->LoadHistograms("Multiplicity");
1163
1164   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1165   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1166   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1167
1168   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1169   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1170
1171   DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1172 }
1173
1174 void SystematicMisalignmentTPC()
1175 {
1176   gSystem->Load("libPWG0base");
1177
1178   SetTPC();
1179
1180   TFile::Open(correctionFile);
1181   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1182   mult->LoadHistograms("Multiplicity");
1183
1184   TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
1185   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1186   mult2->LoadHistograms("Multiplicity");
1187
1188   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1189   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1190   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1191
1192   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1193   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1194
1195   DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
1196 }
1197
1198 void EfficiencySpecies()
1199 {
1200   gSystem->Load("libPWG0base");
1201
1202   Int_t marker[] = {24, 25, 26};
1203   Int_t color[] = {1, 2, 4};
1204
1205   // SPD TPC
1206   const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
1207   Float_t etaRange[] = {0.49, 0.9};
1208   const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
1209
1210   TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
1211   canvas->Divide(2, 1);
1212
1213   for (Int_t loop=0; loop<2; ++loop)
1214   {
1215     Printf("%s", fileName[loop]);
1216
1217     AliCorrection* correction[4];
1218
1219     canvas->cd(loop+1);
1220
1221     gPad->SetGridx();
1222     gPad->SetGridy();
1223     gPad->SetRightMargin(0.05);
1224     //gPad->SetTopMargin(0.05);
1225
1226     TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
1227     legend->SetFillColor(0);
1228     legend->SetEntrySeparation(0.2);
1229
1230     Float_t below = 0;
1231     Float_t total = 0;
1232
1233     TFile* file = TFile::Open(fileName[loop]);
1234     if (!file)
1235     {
1236       Printf("Could not open %s", fileName[loop]);
1237       return;
1238     }
1239
1240     for (Int_t i=0; i<3; ++i)
1241     {
1242       Printf("correction %d", i);
1243
1244       TString name; name.Form("correction_%d", i);
1245       correction[i] = new AliCorrection(name, name);
1246       correction[i]->LoadHistograms();
1247
1248       TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1249       TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1250
1251       // limit vtx axis
1252       gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
1253       meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
1254
1255       // 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)
1256       /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1257         for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1258         {
1259           gene->SetBinContent(x, 0, z, 0);
1260           gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1261           meas->SetBinContent(x, 0, z, 0);
1262           meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1263         }*/
1264
1265       // limit eta axis
1266       gene->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1267       meas->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
1268
1269       TH1* genePt = gene->Project3D(Form("z_%d", i));
1270       TH1* measPt = meas->Project3D(Form("z_%d", i));
1271
1272       genePt->Sumw2();
1273       measPt->Sumw2();
1274
1275       TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1276       effPt->Reset();
1277       effPt->Divide(measPt, genePt, 1, 1, "B");
1278
1279       Int_t bin = 0;
1280       for (bin=20; bin>=1; bin--)
1281       {
1282         if (effPt->GetBinContent(bin) < 0.5)
1283           break;
1284       }
1285
1286       Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1287
1288       Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1289       Printf("%.4f of the particles are below that momentum", fraction);
1290
1291       below += genePt->Integral(1, bin);
1292       total += genePt->Integral();
1293
1294       effPt->SetLineColor(color[i]);
1295       effPt->SetMarkerColor(color[i]);
1296       effPt->SetMarkerStyle(marker[i]);
1297
1298       effPt->GetXaxis()->SetRangeUser(0.06, 1);
1299       effPt->GetYaxis()->SetRangeUser(0, 1);
1300
1301       effPt->GetYaxis()->SetTitleOffset(1.2);
1302
1303       effPt->SetStats(kFALSE);
1304       effPt->SetTitle(titles[loop]);
1305       effPt->GetYaxis()->SetTitle("Efficiency");
1306
1307       effPt->DrawCopy((i == 0) ? "" : "SAME");
1308
1309       legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")));
1310     }
1311
1312     Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
1313
1314     legend->Draw();
1315   }
1316
1317   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1318 }
1319
1320 void ParticleSpeciesComparison1(Bool_t chi2 = kTRUE, const char* fileNameMC = "multiplicityMC_400k_syst_species.root", const char* fileNameESD = "multiplicityMC_100k_syst.root")
1321 {
1322   gSystem->Load("libPWG0base");
1323
1324   TFile::Open(fileNameESD);
1325   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
1326   TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
1327
1328   TH1* results[10];
1329
1330   // loop over cases (normal, enhanced/reduced ratios)
1331   Int_t nMax = 7;
1332   for (Int_t i = 0; i<nMax; ++i)
1333   {
1334     TString folder;
1335     folder.Form("Multiplicity_%d", i);
1336
1337     AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
1338
1339     TFile::Open(fileNameMC);
1340     mult->LoadHistograms();
1341
1342     mult->SetMultiplicityESD(etaRange, hist);
1343
1344     if (chi2)
1345     {
1346       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1347       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1348       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1349     }
1350     else
1351     {
1352       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1353       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1354     }
1355
1356     //Float_t averageRatio = 0;
1357     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1358
1359     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1360
1361     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1362   }
1363
1364   DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
1365
1366   TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
1367
1368   for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1369   {
1370     results[0]->SetBinError(i, 0);
1371     mc->SetBinError(i, 0);
1372   }
1373
1374   const char* legendStrings[] = { "#pi^{#pm}", 0, "K^{#pm}", 0, "p,#bar{p}", 0 };
1375
1376   DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
1377
1378   //not valid: draw chi2 uncertainty on top!
1379   /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1380   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1381   errorHist->SetLineColor(1);
1382   errorHist->SetLineWidth(2);
1383   TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1384   for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1385   {
1386     errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1387     errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1388   }
1389
1390   errorHist->DrawCopy("SAME");
1391   errorHist2->DrawCopy("SAME");*/
1392
1393   //canvas->SaveAs(canvas->GetName());
1394
1395   DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
1396
1397   //errorHist->DrawCopy("SAME");
1398   //errorHist2->DrawCopy("SAME");
1399
1400   //canvas2->SaveAs(canvas2->GetName());
1401 }
1402
1403 /*void ParticleSpeciesComparison2()
1404 {
1405   gSystem->Load("libPWG0base");
1406
1407   const char* fileNameMC = "multiplicityMC_400k_syst.root";
1408   const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1409   Bool_t chi2 = 0;
1410
1411   TFile::Open(fileNameMC);
1412   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1413   mult->LoadHistograms();
1414
1415   TH1* mc[10];
1416   TH1* results[10];
1417
1418   // loop over cases (normal, enhanced/reduced ratios)
1419   Int_t nMax = 7;
1420   for (Int_t i = 0; i<nMax; ++i)
1421   {
1422     TString folder;
1423     folder.Form("Multiplicity_%d", i);
1424
1425     AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1426
1427     TFile::Open(fileNameESD);
1428     mult2->LoadHistograms();
1429
1430     mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1431
1432     if (chi2)
1433     {
1434       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1435       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1436       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1437     }
1438     else
1439     {
1440       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1441       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1442     }
1443
1444     //Float_t averageRatio = 0;
1445     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1446
1447     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1448
1449     TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1450     mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1451
1452     //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1453     //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1454
1455     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1456   }
1457
1458   DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1459 }*/
1460
1461 TH1* Invert(TH1* eff)
1462 {
1463   // calculate corr = 1 / eff
1464
1465   TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1466   corr->Reset();
1467
1468   for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1469   {
1470     if (eff->GetBinContent(i) > 0)
1471     {
1472       corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1473       corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1474     }
1475   }
1476
1477   return corr;
1478 }
1479
1480 void TriggerVertexCorrection()
1481 {
1482   //
1483   // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1484   //
1485
1486   gSystem->Load("libPWG0base");
1487
1488   TFile::Open(correctionFile);
1489   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1490   mult->LoadHistograms("Multiplicity");
1491
1492   TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1493   TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1494
1495   TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
1496
1497   corrINEL->SetStats(kFALSE);
1498   corrINEL->GetXaxis()->SetRangeUser(0, 30);
1499   corrINEL->GetYaxis()->SetRangeUser(0, 5);
1500   corrINEL->SetTitle(";true multiplicity;correction factor");
1501   corrINEL->SetMarkerStyle(22);
1502   corrINEL->Draw("PE");
1503
1504   corrMB->SetStats(kFALSE);
1505   corrMB->SetLineColor(2);
1506   corrMB->SetMarkerStyle(25);
1507   corrMB->SetMarkerColor(2);
1508   corrMB->Draw("SAME PE");
1509
1510   TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
1511   legend->SetFillColor(0);
1512   legend->AddEntry(corrINEL, "correction to inelastic sample");
1513   legend->AddEntry(corrMB, "correction to minimum bias sample");
1514
1515   legend->Draw();
1516
1517   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1518 }
1519
1520 void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
1521 {
1522   gSystem->Load("libPWG0base");
1523
1524   TFile::Open(correctionFile);
1525   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1526   mult->LoadHistograms("Multiplicity");
1527
1528   TFile::Open(measuredFile);
1529   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1530   mult2->LoadHistograms("Multiplicity");
1531
1532   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1533
1534   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1535
1536   TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
1537
1538   TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1539   TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
1540
1541   if (!mc)
1542   {
1543     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1544     DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
1545   }
1546
1547   TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
1548   canvas->SetGridx();
1549   canvas->SetGridy();
1550   canvas->SetRightMargin(0.05);
1551   canvas->SetTopMargin(0.05);
1552
1553   errorResponse->SetLineColor(1);
1554   errorResponse->GetXaxis()->SetRangeUser(0, 200);
1555   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1556   errorResponse->SetStats(kFALSE);
1557   errorResponse->SetTitle(";true multiplicity;Uncertainty");
1558
1559   errorResponse->Draw();
1560
1561   errorMeasured->SetLineColor(2);
1562   errorMeasured->Draw("SAME");
1563
1564   errorBoth->SetLineColor(4);
1565   errorBoth->Draw("SAME");
1566
1567   Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
1568   Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
1569   Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
1570
1571   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1572
1573   TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
1574   errorResponse->Write();
1575   errorMeasured->Write();
1576   errorBoth->Write();
1577   file->Close();
1578 }
1579
1580 void StatisticalUncertaintyCompare(const char* det = "SPD")
1581 {
1582   TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
1583   TH1* errorResponse = (TH1*) file1->Get("errorResponse");
1584   TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
1585   TH1* errorBoth = (TH1*) file1->Get("errorBoth");
1586
1587   TString str;
1588   str.Form("StatisticalUncertaintyCompare%s", det);
1589
1590   TCanvas* canvas = new TCanvas(str, str, 600, 400);
1591   canvas->SetGridx();
1592   canvas->SetGridy();
1593   canvas->SetRightMargin(0.05);
1594   canvas->SetTopMargin(0.05);
1595
1596   errorResponse->SetLineColor(1);
1597   errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
1598   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1599   errorResponse->SetStats(kFALSE);
1600   errorResponse->SetTitle(";true multiplicity;Uncertainty");
1601
1602   errorResponse->Draw();
1603
1604   errorMeasured->SetLineColor(2);
1605   errorMeasured->Draw("SAME");
1606
1607   errorBoth->SetLineColor(4);
1608   errorBoth->Draw("SAME");
1609
1610   TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
1611   TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
1612
1613   errorBoth2->SetLineColor(4);
1614   errorBoth2->SetLineStyle(2);
1615   errorBoth2->Draw("SAME");
1616
1617   TLegend* legend = new TLegend(0.2, 0.6, 0.6, 0.9);
1618   legend->SetFillColor(0);
1619   legend->AddEntry(errorResponse, "response matrix (Bayesian)");
1620   legend->AddEntry(errorMeasured, "measured (Bayesian)");
1621   legend->AddEntry(errorBoth, "both (Bayesian)");
1622   legend->AddEntry(errorBoth2, "both (#chi^{2} minimization)");
1623   legend->Draw();
1624
1625   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1626 }
1627
1628 void EfficiencyComparison(Int_t eventType = 2)
1629 {
1630  const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
1631
1632   gSystem->Load("libPWG0base");
1633
1634   TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
1635   canvas->SetGridx();
1636   canvas->SetGridy();
1637   canvas->SetRightMargin(0.05);
1638   canvas->SetTopMargin(0.05);
1639
1640   AliMultiplicityCorrection* data[4];
1641   TH1* effArray[4];
1642
1643   Int_t markers[] = { 24, 25, 26, 5 };
1644   Int_t colors[] = { 1, 2, 3, 4 };
1645
1646   TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
1647   legend->SetFillColor(0);
1648
1649   TH1* effError = 0;
1650
1651   for (Int_t i=0; i<4; ++i)
1652   {
1653     TString name;
1654     name.Form("Multiplicity_%d", i);
1655
1656     TFile::Open(files[i]);
1657     data[i] = new AliMultiplicityCorrection(name, name);
1658
1659     if (i < 3)
1660     {
1661       data[i]->LoadHistograms("Multiplicity");
1662     }
1663     else
1664       data[i]->LoadHistograms("Multiplicity_0");
1665
1666     TH1* eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
1667     effArray[i] = eff;
1668
1669     eff->GetXaxis()->SetRangeUser(0, 15);
1670     eff->GetYaxis()->SetRangeUser(0, 1.1);
1671     eff->SetStats(kFALSE);
1672     eff->SetTitle(";true multiplicity;Efficiency");
1673     eff->SetLineColor(colors[i]);
1674     eff->SetMarkerColor(colors[i]);
1675     eff->SetMarkerStyle(markers[i]);
1676
1677     if (i == 3)
1678     {
1679       for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1680         eff->SetBinError(bin, 0);
1681
1682       // loop over cross section combinations
1683       for (Int_t j=1; j<7; ++j)
1684       {
1685         AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
1686         mult->LoadHistograms(Form("Multiplicity_%d", j));
1687
1688         TH1* eff2 = mult->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType);
1689
1690         for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1691         {
1692           // TODO we could also do asymmetric errors here
1693           Float_t deviation = TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
1694
1695           eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), (Double_t) deviation));
1696         }
1697       }
1698
1699       for (Int_t bin=1; bin<=20; bin++)
1700         if (eff->GetBinContent(bin) > 0)
1701           Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1702
1703       effError = (TH1*) eff->Clone("effError");
1704       effError->Reset();
1705
1706       for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
1707         if (eff->GetBinContent(bin) > 0)
1708           effError->SetBinContent(bin, 10.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1709
1710       effError->SetLineColor(1);
1711       effError->SetMarkerStyle(1);
1712       effError->DrawCopy("SAME HIST");
1713     }
1714
1715     eff->SetBinContent(1, 0);
1716     eff->SetBinError(1, 0);
1717
1718     canvas->cd();
1719     if (i == 0)
1720     {
1721       eff->DrawCopy("P");
1722     }
1723     else
1724       eff->DrawCopy("SAME P");
1725
1726     legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
1727   }
1728
1729   legend->AddEntry(effError, "relative syst. uncertainty #times 10");
1730
1731   legend->Draw();
1732
1733   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1734 }
1735
1736 void ModelDependencyPlot()
1737 {
1738   gSystem->Load("libPWG0base");
1739
1740   TFile::Open("multiplicityMC_3M.root");
1741   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1742   mult->LoadHistograms("Multiplicity");
1743
1744   TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
1745
1746   TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
1747   canvas->SetGridx();
1748   canvas->SetGridy();
1749   //canvas->SetRightMargin(0.05);
1750   //canvas->SetTopMargin(0.05);
1751
1752   canvas->Divide(2, 1);
1753
1754   canvas->cd(2);
1755   gPad->SetLogy();
1756
1757   Int_t selectedMult = 30;
1758   Int_t yMax = 200000;
1759
1760   TH1* full = proj->ProjectionX("full");
1761   TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
1762
1763   full->SetStats(kFALSE);
1764   full->GetXaxis()->SetRangeUser(0, 200);
1765   full->GetYaxis()->SetRangeUser(5, yMax);
1766   full->SetTitle(";multiplicity");
1767
1768   selected->SetLineColor(0);
1769   selected->SetMarkerColor(2);
1770   selected->SetMarkerStyle(7);
1771
1772   full->Draw();
1773   selected->Draw("SAME P");
1774
1775   TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
1776   legend->SetFillColor(0);
1777   legend->AddEntry(full, "true");
1778   legend->AddEntry(selected, "measured");
1779   legend->Draw();
1780
1781   TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
1782   line->SetLineWidth(2);
1783   line->Draw();
1784
1785   canvas->cd(1);
1786   gPad->SetLogy();
1787
1788   full = proj->ProjectionY("full2");
1789   selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
1790
1791   full->SetStats(kFALSE);
1792   full->GetXaxis()->SetRangeUser(0, 200);
1793   full->GetYaxis()->SetRangeUser(5, yMax);
1794   full->SetTitle(";multiplicity");
1795
1796   full->SetLineColor(0);
1797   full->SetMarkerColor(2);
1798   full->SetMarkerStyle(7);
1799
1800   full->Draw("P");
1801   selected->Draw("SAME");
1802
1803   legend->Draw();
1804
1805   line = new TLine(selectedMult, 5, selectedMult, yMax);
1806   line->SetLineWidth(2);
1807   line->Draw();
1808
1809   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1810 }
1811
1812 void SystematicpTSpectrum()
1813 {
1814   gSystem->Load("libPWG0base");
1815
1816   TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
1817   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1818   mult->LoadHistograms("Multiplicity");
1819
1820   TFile::Open("multiplicityMC_100k_syst.root");
1821   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1822   mult2->LoadHistograms("Multiplicity");
1823
1824   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1825   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1826   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1827
1828   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1829   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1830
1831   DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
1832 }
1833
1834 // to be deleted
1835 /*void covMatrix(Bool_t mc = kTRUE)
1836 {
1837   gSystem->Load("libPWG0base");
1838
1839   TFile::Open(correctionFile);
1840   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1841   mult->LoadHistograms("Multiplicity");
1842
1843   TFile::Open(measuredFile);
1844   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1845   mult2->LoadHistograms("Multiplicity");
1846
1847   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1848
1849   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1850
1851   mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
1852 }*/
1853
1854 Double_t FitPtFunc(Double_t *x, Double_t *par)
1855 {
1856   Double_t xx = x[0];
1857
1858   Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
1859   Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
1860
1861   const Float_t kTransitionWidth = 0;
1862
1863   // power law part
1864   if (xx < par[0] - kTransitionWidth)
1865   {
1866     return val1;
1867   }
1868   /*else if (xx < par[0] + kTransitionWidth)
1869   {
1870     // smooth transition
1871     Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
1872     return (1 - factor) * val1 + factor * val2;
1873   }*/
1874   else
1875   {
1876     return val2;
1877   }
1878 }
1879
1880 void FitPt(const char* fileName = "firstplots100k_truept.root")
1881 {
1882   gSystem->Load("libPWG0base");
1883
1884   TFile::Open(fileName);
1885
1886   /*
1887   // merge corrections
1888   AliCorrection* correction[4];
1889   TList list;
1890
1891   for (Int_t i=0; i<4; ++i)
1892   {
1893     Printf("correction %d", i);
1894
1895     TString name; name.Form("correction_%d", i);
1896     correction[i] = new AliCorrection(name, name);
1897     correction[i]->LoadHistograms();
1898
1899     if (i > 0)
1900       list.Add(correction[i]);
1901   }
1902
1903   correction[0]->Merge(&list);
1904
1905   TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
1906
1907   // limit vtx, eta axis
1908   gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
1909   gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
1910
1911   TH1* genePt = gene->Project3D("z");*/
1912   TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
1913   genePt->Sumw2();
1914
1915   //genePt->Scale(1.0 / genePt->Integral());
1916
1917   // normalize by bin width
1918   for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1919     genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1920
1921   /// genePt->GetXaxis()->GetBinCenter(x));
1922
1923   genePt->GetXaxis()->SetRangeUser(0, 7.9);
1924   //genePt->GetYaxis()->SetTitle("a.u.");
1925
1926   TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
1927   //func->SetLineColor(2);
1928   func->SetParameters(1, -1, 1, 1, 1);
1929   func->SetParLimits(3, 1, 10);
1930   func->SetParLimits(4, 0, 10);
1931
1932   //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
1933
1934   //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
1935   //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
1936   //func->FixParameter(0, 0.314);
1937   //func->SetParLimits(0, 0.1, 0.3);
1938
1939   genePt->SetMarkerStyle(25);
1940   genePt->SetTitle("");
1941   genePt->SetStats(kFALSE);
1942   genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
1943   //func->Draw("SAME");
1944
1945   // fit only exp. part
1946   func->SetParameters(1, -1);
1947   func->FixParameter(2, 0);
1948   func->FixParameter(3, 1);
1949   func->FixParameter(4, 1);
1950   genePt->Fit(func, "0", "", 0.2, 1);
1951
1952   new TCanvas;
1953   genePt->DrawCopy("P");
1954   func->SetRange(0.02, 8);
1955   func->DrawCopy("SAME");
1956   gPad->SetLogy();
1957
1958   // now fix exp. parameters and fit second part
1959   Double_t param0 = func->GetParameter(0);
1960   Double_t param1 = func->GetParameter(1);
1961   func->SetParameters(0, -1, 1, 1, 1);
1962   func->FixParameter(0, 0);
1963   func->FixParameter(1, -1);
1964   func->ReleaseParameter(2);
1965   func->ReleaseParameter(3);
1966   func->ReleaseParameter(4);
1967   func->SetParLimits(3, 1, 10);
1968   func->SetParLimits(4, 0, 10);
1969
1970   genePt->Fit(func, "0", "", 1.5, 4);
1971
1972   new TCanvas;
1973   genePt->DrawCopy("P");
1974   func->SetRange(0.02, 8);
1975   func->DrawCopy("SAME");
1976   gPad->SetLogy();
1977
1978   // fit both
1979   func->SetParameter(0, param0);
1980   func->SetParameter(1, param1);
1981   func->ReleaseParameter(0);
1982   func->ReleaseParameter(1);
1983
1984   new TCanvas;
1985   genePt->DrawCopy("P");
1986   func->SetRange(0.02, 5);
1987   func->DrawCopy("SAME");
1988   gPad->SetLogy();
1989
1990   genePt->Fit(func, "0", "", 0.2, 4);
1991
1992   TCanvas* canvas = new TCanvas("FitPt", "FitPt", 800, 400);
1993   canvas->Divide(2, 1);
1994   canvas->cd(1);
1995
1996   gPad->SetGridx();
1997   gPad->SetGridy();
1998   gPad->SetLeftMargin(0.13);
1999   gPad->SetRightMargin(0.05);
2000   gPad->SetTopMargin(0.05);
2001
2002   genePt->GetXaxis()->SetRangeUser(0, 4.9);
2003   genePt->GetYaxis()->SetRangeUser(1e-2, 1e4);
2004   genePt->GetYaxis()->SetTitleOffset(1.4);
2005   genePt->GetXaxis()->SetTitleOffset(1.1);
2006   genePt->DrawCopy("P");
2007   func->SetRange(0.02, 5);
2008   func->DrawCopy("SAME");
2009   gPad->SetLogy();
2010
2011   canvas->cd(2);
2012
2013   TH1* genePtClone = (TH1*) genePt->Clone("genePtClone");
2014   genePtClone->Reset();
2015   genePtClone->DrawCopy("P");
2016
2017   gPad->SetGridx();
2018   gPad->SetGridy();
2019   gPad->SetLeftMargin(0.13);
2020   gPad->SetRightMargin(0.05);
2021   gPad->SetTopMargin(0.05);
2022
2023   func->DrawCopy("SAME");
2024   gPad->SetLogy();
2025
2026   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2027
2028   TH1* first = (TH1*) func->GetHistogram()->Clone("first");
2029
2030   TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
2031
2032   TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
2033
2034   for (Int_t param=0; param<5; param++)
2035   {
2036     for (Int_t sign=0; sign<2; sign++)
2037     {
2038       TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign));  // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
2039       func2->SetParameters(func->GetParameters());
2040       //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
2041
2042       Float_t factor = ((sign == 0) ? 0.9 : 1.1);
2043       func2->SetParameter(param, func2->GetParameter(param) * factor);
2044       //func2->Print();
2045
2046       canvas->cd(2);
2047       func2->SetLineWidth(1);
2048       func2->SetLineColor(2);
2049       func2->DrawCopy("SAME");
2050
2051       canvas2->cd();
2052       TH1* second = func2->GetHistogram();
2053       second->Divide(first);
2054       second->SetLineColor(param + 1);
2055       second->GetYaxis()->SetRangeUser(0, 2);
2056       second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
2057       second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
2058     }
2059   }
2060
2061   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
2062   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
2063
2064   file->Close();
2065 }
2066
2067 void DrawSystematicpT()
2068 {
2069   TFile* file = TFile::Open("SystematicpT.root");
2070
2071   TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
2072   TH1* result2 = (TH1*) file->Get("result_unity");
2073
2074   TH1* mcHist[12];
2075   TH1* result[12];
2076
2077   Int_t nParams = 5;
2078
2079   for (Int_t id=0; id<nParams*2; ++id)
2080   {
2081     mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
2082     result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
2083   }
2084
2085   DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
2086
2087   //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2088
2089   DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
2090
2091   //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
2092
2093   // does not make sense: mc is different
2094   //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
2095 }
2096
2097 void SystematicpT(Bool_t chi2 = 1)
2098 {
2099   gSystem->Load("libPWG0base");
2100
2101   TFile::Open("ptspectrum900.root");
2102   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2103   mult->LoadHistograms("Multiplicity");
2104
2105   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2106
2107   TH1* mcHist[12];
2108   TH1* result[12];
2109
2110   Int_t nParams = 5;
2111
2112   for (Int_t param=0; param<nParams; param++)
2113   {
2114     for (Int_t sign=0; sign<2; sign++)
2115     {
2116       // calculate result with systematic effect
2117       TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
2118       mult2->LoadHistograms("Multiplicity");
2119
2120       mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2121
2122       if (chi2)
2123       {
2124         mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2125         mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2126       }
2127       else
2128         mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2129
2130       Int_t id = param * 2 + sign;
2131
2132       mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
2133       result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
2134
2135       TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
2136       DrawResultRatio(mcHist[id], result[id], tmp);
2137     }
2138   }
2139
2140   // calculate normal result
2141   TFile::Open("ptspectrum100_1.root");
2142   mult2->LoadHistograms("Multiplicity");
2143   TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
2144
2145   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2146
2147   if (chi2)
2148   {
2149     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2150   }
2151   else
2152     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2153
2154   TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
2155
2156   TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
2157   mcHist2->Write();
2158   result2->Write();
2159   for (Int_t id=0; id<nParams*2; ++id)
2160   {
2161     mcHist[id]->Write();
2162     result[id]->Write();
2163   }
2164   file->Close();
2165
2166   DrawSystematicpT();
2167 }
2168
2169 void DrawSystematicpT2()
2170 {
2171   //displayRange = 200;
2172
2173   // read from file
2174   TFile* file = TFile::Open("SystematicpT2.root");
2175   TH1* mcHist = (TH1*) file->Get("mymc");
2176   TH1* result[12];
2177   result[0] = (TH1*) file->Get("result_unity");
2178   Int_t nParams = 5;
2179   for (Int_t id=0; id<nParams*2; ++id)
2180     result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
2181
2182   DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
2183   DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
2184   DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
2185 }
2186
2187 void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2188 {
2189   gSystem->Load("libPWG0base");
2190
2191   if (tpc)
2192   {
2193     SetTPC();
2194     TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2195   }
2196   else
2197     TFile::Open("ptspectrum100_1.root");
2198
2199   AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2200   measured->LoadHistograms("Multiplicity");
2201   TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2202
2203   TH1* result[12];
2204
2205   Int_t nParams = 5;
2206
2207   // -1 = unity change, 0...4 parameters
2208   for (Int_t id=-1; id<nParams*2; id++)
2209   {
2210     Int_t param = id / 2;
2211     Int_t sign = id % 2;
2212
2213     TString idStr;
2214     if (id == -1)
2215     {
2216       idStr = "unity";
2217     }
2218     else
2219       idStr.Form("%d_%d", param, sign);
2220
2221     // calculate result with systematic effect
2222     if (tpc)
2223     {
2224       TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
2225     }
2226     else
2227       TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
2228
2229     AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2230     response->LoadHistograms("Multiplicity");
2231
2232     response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
2233
2234     if (chi2)
2235     {
2236       response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2237       response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2238     }
2239     else
2240       response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
2241
2242     result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
2243
2244     TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
2245     DrawResultRatio(mcHist, result[id+1], tmp);
2246   }
2247
2248   TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
2249   mcHist->Write();
2250   for (Int_t id=0; id<nParams*2+1; ++id)
2251     result[id]->Write();
2252   file->Close();
2253
2254   DrawSystematicpT2();
2255 }
2256
2257 void SystematicpTCutOff(Bool_t chi2 = kTRUE)
2258 {
2259   // only needed for TPC
2260   SetTPC();
2261
2262   gSystem->Load("libPWG0base");
2263
2264   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
2265   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2266   mult->LoadHistograms("Multiplicity");
2267
2268   TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
2269   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2270   mult2->LoadHistograms("Multiplicity");
2271
2272   // "normal" result
2273   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2274
2275   if (chi2)
2276   {
2277     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2278     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2279   }
2280   else
2281     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2282
2283   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
2284   TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
2285
2286   TH1* syst[2];
2287
2288   // change of pt spectrum (down)
2289   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
2290   AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
2291   mult3->LoadHistograms("Multiplicity");
2292   mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2293   if (chi2)
2294   {
2295     mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2296     mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2297   }
2298   else
2299     mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2300   syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
2301
2302   // change of pt spectrum (up)
2303   TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
2304   AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
2305   mult4->LoadHistograms("Multiplicity");
2306   mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2307   if (chi2)
2308   {
2309     mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2310     mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
2311   }
2312   else
2313     mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
2314   syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
2315
2316   DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
2317
2318   Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
2319   Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
2320 }
2321
2322 TH1* SystematicsSummary(Bool_t tpc = 1)
2323 {
2324   Int_t nEffects = 7;
2325
2326   TH1* effects[10];
2327   const char** names = 0;
2328   Int_t colors[] = { 1, 2, 3, 4, 6, 7, 8 };
2329   Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
2330
2331   for (Int_t i=0; i<nEffects; ++i)
2332     effects[i] = new TH1F("SystematicsSummary", ";true multiplicity;Effect", 201, -0.5, 200.5);
2333
2334   if (tpc)
2335   {
2336     SetTPC();
2337
2338     const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
2339     names = namesTPC;
2340
2341     // method
2342     TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
2343     TH1* hist = (TH1*) file->Get("errorBoth");
2344
2345     // smooth a bit, but skip 0 bin
2346     effects[0]->SetBinContent(2, hist->GetBinContent(2));
2347     for (Int_t i=3; i<=200; ++i)
2348       effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2349
2350     // relative x-section
2351     effects[1]->SetBinContent(2, 0.005);
2352     effects[1]->SetBinContent(3, 0.0025);
2353     effects[1]->SetBinContent(4, 0.0025);
2354
2355     // particle composition
2356     for (Int_t i=2; i<=101; ++i)
2357     {
2358       if (i < 41)
2359       {
2360         effects[2]->SetBinContent(i, 0.01);
2361       }
2362       else if (i < 76)
2363       {
2364         effects[2]->SetBinContent(i, 0.02);
2365       }
2366       else
2367         effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
2368     }
2369
2370     // pt cut off (only tpc)
2371     for (Int_t i=2; i<=101; ++i)
2372     {
2373       if (i < 11)
2374       {
2375         effects[3]->SetBinContent(i, 0.05);
2376       }
2377       else if (i < 51)
2378       {
2379         effects[3]->SetBinContent(i, 0.01);
2380       }
2381       else
2382         effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
2383     }
2384
2385     // track selection (only tpc)
2386     for (Int_t i=2; i<=101; ++i)
2387       effects[4]->SetBinContent(i, 0.03);
2388
2389     // secondaries
2390     for (Int_t i=2; i<=101; ++i)
2391       effects[5]->SetBinContent(i, 0.01);
2392
2393     // pt spectrum
2394     for (Int_t i=2; i<=101; ++i)
2395     {
2396       if (i < 21)
2397       {
2398         effects[6]->SetBinContent(i, 0.05);
2399       }
2400       else if (i < 51)
2401       {
2402         effects[6]->SetBinContent(i, 0.02);
2403       }
2404       else
2405         effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
2406     }
2407
2408   }
2409   else
2410   {
2411     displayRange = 200;
2412     nEffects = 5;
2413
2414     const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Rel. cross-section", "Particle composition", "Secondaries", "p_{t} spectrum"};
2415     names = namesSPD;
2416
2417     // method
2418     TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
2419     TH1* hist = (TH1*) file->Get("errorBoth");
2420
2421     // smooth a bit, but skip 0 bin
2422     effects[0]->SetBinContent(2, hist->GetBinContent(2));
2423     for (Int_t i=3; i<=201; ++i)
2424       effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
2425
2426     // relative x-section
2427     effects[1]->SetBinContent(2, 0.01);
2428     effects[1]->SetBinContent(3, 0.005);
2429
2430     // particle composition
2431     for (Int_t i=2; i<=201; ++i)
2432     {
2433       if (i < 6)
2434       {
2435         effects[2]->SetBinContent(i, 0.3);
2436       }
2437       else if (i < 11)
2438       {
2439         effects[2]->SetBinContent(i, 0.05);
2440       }
2441       else if (i < 121)
2442       {
2443         effects[2]->SetBinContent(i, 0.02);
2444       }
2445       else if (i < 151)
2446       {
2447         effects[2]->SetBinContent(i, 0.02 + 0.04 / 30 * (i - 121));
2448       }
2449       else
2450         effects[2]->SetBinContent(i, 0.06 + 0.1 / 30 * (i - 151));
2451     }
2452
2453     // secondaries
2454     for (Int_t i=2; i<=201; ++i)
2455       effects[3]->SetBinContent(i, 0.01);
2456
2457     // pt spectrum
2458     for (Int_t i=2; i<=201; ++i)
2459     {
2460       if (i < 6)
2461       {
2462         effects[4]->SetBinContent(i, 1);
2463       }
2464       else if (i < 121)
2465       {
2466         effects[4]->SetBinContent(i, 0.03);
2467       }
2468       else if (i < 151)
2469       {
2470         effects[4]->SetBinContent(i, 0.03 + 0.07 / 30 * (i - 121));
2471       }
2472       else
2473         effects[4]->SetBinContent(i, 0.1);
2474     }
2475   }
2476
2477   TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 400);
2478   canvas->SetRightMargin(0.25);
2479   canvas->SetTopMargin(0.05);
2480   TLegend* legend = new TLegend(0.2, 0.4, 0.5, 0.4 + 0.5 * nEffects / 7);
2481   legend->SetFillColor(0);
2482
2483   for (Int_t i=0; i<nEffects; ++i)
2484   {
2485     TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
2486     /*current->Reset();
2487     for (Int_t j=0; j<nEffects-i; ++j)
2488       current->Add(effects[j]);*/
2489
2490     current->SetLineColor(colors[i]);
2491     //current->SetFillColor(colors[i]);
2492     current->SetMarkerColor(colors[i]);
2493     //current->SetMarkerStyle(markers[i]);
2494
2495     current->SetStats(kFALSE);
2496     current->GetYaxis()->SetRangeUser(0, 0.4);
2497     current->GetXaxis()->SetRangeUser(0, displayRange);
2498     current->DrawCopy(((i == 0) ? "" : "SAME"));
2499     legend->AddEntry(current, names[i]);
2500
2501     TLatex* text = new TLatex(displayRange+5, current->GetBinContent(displayRange+1), names[i]);
2502     text->SetTextColor(colors[i]);
2503     text->Draw();
2504   }
2505
2506   // add total in square
2507   TH1* total = (TH1*) effects[0]->Clone("total");
2508   total->Reset();
2509
2510   for (Int_t i=0; i<nEffects; ++i)
2511   {
2512     //Printf("%d %f", i, effects[i]->GetBinContent(20));
2513     effects[i]->Multiply(effects[i]);
2514     total->Add(effects[i]);
2515   }
2516
2517   for (Int_t i=1; i<=total->GetNbinsX(); ++i)
2518     if (total->GetBinContent(i) > 0)
2519       total->SetBinContent(i, TMath::Min(sqrt(total->GetBinContent(i)), 1.0));
2520
2521   //Printf("%f", total->GetBinContent(20));
2522
2523   total->SetMarkerStyle(3);
2524   total->SetMarkerColor(1);
2525   legend->AddEntry(total, "total");
2526   total->DrawCopy("SAME P");
2527
2528   legend->Draw();
2529
2530   canvas->SaveAs(canvas->GetName());
2531
2532   return total;
2533 }
2534
2535 void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
2536 {
2537   gSystem->Load("libPWG0base");
2538
2539   if (tpc)
2540     SetTPC();
2541
2542   if (!chi2)
2543     Printf("WARNING: Bayesian set. This is only for test!");
2544
2545   // systematic error
2546   TH1* error = SystematicsSummary(tpc);
2547
2548   TFile::Open(correctionFile);
2549   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2550   mult->LoadHistograms("Multiplicity");
2551
2552   TFile::Open(measuredFile);
2553   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
2554   mult2->LoadHistograms("Multiplicity");
2555
2556   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
2557
2558   if (chi2)
2559   {
2560     mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
2561     mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kINEL);
2562   }
2563   else
2564     mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kINEL, 1, 100, 0, kFALSE);
2565
2566   TH1* mcHist = mult2->GetMultiplicityINEL(etaRange)->ProjectionY("mymc");
2567   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
2568
2569   DrawResultRatio(mcHist, result, "finalPlotCheck.eps");
2570
2571   // normalize result
2572   result->Scale(1.0 / result->Integral(2, 200));
2573
2574   result->GetXaxis()->SetRangeUser(0, ((tpc) ? displayRange : 200));
2575   result->SetBinContent(1, 0); result->SetBinError(1, 0);
2576   result->SetTitle(";true multiplicity;Probability");
2577   result->SetLineColor(1);
2578   result->SetStats(kFALSE);
2579
2580   TH1* systError = (TH1*) result->Clone("systError");
2581   for (Int_t i=2; i<=systError->GetNbinsX(); ++i)
2582     systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
2583
2584   // change error drawing style
2585   systError->SetFillColor(15);
2586
2587   TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", 800, 400);
2588   canvas->SetRightMargin(0.05);
2589   canvas->SetTopMargin(0.05);
2590
2591   systError->Draw("E2 ][");
2592   result->DrawCopy("SAME E ][");
2593   canvas->SetLogy();
2594
2595   //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
2596   TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
2597   text->SetFillColor(0);
2598   text->SetTextAlign(12);
2599   text->AddText("Systematic errors summed quadratically");
2600   text->AddText("0.6 million minimum bias events");
2601   text->AddText("corrected to inelastic events");
2602   text->Draw("B");
2603
2604   TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
2605   text2->SetFillColor(0);
2606   text2->SetTextAlign(12);
2607   text2->AddText("#sqrt{s} = 14 TeV");
2608   if (tpc)
2609   {
2610     text2->AddText("|#eta| < 0.9");
2611   }
2612   else
2613     text2->AddText("|#eta| < 2.0");
2614   text2->AddText("simulated data (PYTHIA)");
2615   text2->Draw("B");
2616
2617   if (tpc)
2618   {
2619     TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
2620     text3->SetNDC();
2621     text3->Draw();
2622   }
2623   else
2624   {
2625     TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
2626     text3->SetNDC();
2627     text3->Draw();
2628   }
2629
2630   // alice logo
2631   TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
2632   pad->Draw();
2633   pad->cd();
2634   TImage* img = TImage::Open("$HOME/alice.png");
2635   img->SetImageQuality(TAttImage::kImgBest);
2636   img->Draw();
2637
2638   canvas->Modified();
2639
2640 /*  TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
2641   text->SetTextSize(0.04);
2642   text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
2643   text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
2644   text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
2645   text->Draw();*/
2646
2647
2648   canvas->SaveAs(canvas->GetName());
2649 }