PWG0/dNdEta/plotsMultiplicity.C

   1 //
   2 // plots for the note about multplicity measurements
   3 //
   4
   5 const char* correctionFile = "multiplicityMC_2M.root";
   6 const char* measuredFile   = "multiplicityMC_1M_3.root";
   7 Int_t etaRange = 3;
   8
   9 const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
  10 const char* measuredFileTPC   = "multiplicityMC_TPC_0.6M.root";
  11 Int_t etaRangeTPC = 1;
  12
  13 void SetTPC()
  14 {
  15   correctionFile = correctionFileTPC;
  16   measuredFile = measuredFileTPC;
  17   etaRange = etaRangeTPC;
  18 }
  19
  20 void responseMatrixPlot()
  21 {
  22   gSystem->Load("libPWG0base");
  23
  24   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
  25
  26   TFile::Open(correctionFile);
  27   mult->LoadHistograms("Multiplicity");
  28
  29   TH1* hist = mult->GetCorrelation(etaRange)->Project3D("zy");
  30   hist->SetStats(kFALSE);
  31
  32   hist->SetTitle(";true multiplicity;measured multiplicity;Entries");
  33   hist->GetXaxis()->SetRangeUser(0, 200);
  34   hist->GetYaxis()->SetRangeUser(0, 200);
  35
  36   TCanvas* canvas = new TCanvas("c1", "c1", 800, 600);
  37   canvas->SetRightMargin(0.15);
  38   canvas->SetTopMargin(0.05);
  39
  40   gPad->SetLogz();
  41   hist->Draw("COLZ");
  42
  43   canvas->SaveAs("responsematrix.eps");
  44 }
  45
  46 TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
  47 {
  48   // normalize unfolded result to mc hist
  49   result->Scale(1.0 / result->Integral(2, 200));
  50   result->Scale(mcHist->Integral(2, 200));
  51
  52   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
  53   canvas->Range(0, 0, 1, 1);
  54
  55   TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
  56   pad1->Draw();
  57
  58   TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
  59   pad2->Draw();
  60
  61   pad1->SetRightMargin(0.05);
  62   pad2->SetRightMargin(0.05);
  63
  64   // no border between them
  65   pad1->SetBottomMargin(0);
  66   pad2->SetTopMargin(0);
  67
  68   pad1->cd();
  69
  70   mcHist->GetXaxis()->SetLabelSize(0.06);
  71   mcHist->GetYaxis()->SetLabelSize(0.06);
  72   mcHist->GetXaxis()->SetTitleSize(0.06);
  73   mcHist->GetYaxis()->SetTitleSize(0.06);
  74   mcHist->GetYaxis()->SetTitleOffset(0.6);
  75
  76   mcHist->GetXaxis()->SetRangeUser(0, 200);
  77
  78   mcHist->SetTitle(";true multiplicity;Entries");
  79   mcHist->SetStats(kFALSE);
  80
  81   mcHist->DrawCopy("HIST E");
  82   gPad->SetLogy();
  83
  84   result->SetLineColor(2);
  85   result->DrawCopy("SAME HISTE");
  86
  87   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
  88   legend->AddEntry(mcHist, "true distribution");
  89   legend->AddEntry(result, "unfolded distribution");
  90   legend->SetFillColor(0);
  91   legend->Draw();
  92
  93   pad2->cd();
  94   pad2->SetBottomMargin(0.15);
  95
  96   // calculate ratio
  97   mcHist->Sumw2();
  98   TH1* ratio = (TH1*) mcHist->Clone("ratio");
  99   result->Sumw2();
 100   ratio->Divide(ratio, result, 1, 1, "");
 101   ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
 102   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
 103
 104   ratio->DrawCopy();
 105
 106   // get average of ratio
 107   Float_t sum = 0;
 108   for (Int_t i=2; i<=150; ++i)
 109   {
 110     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
 111   }
 112   sum /= 149;
 113
 114   printf("Average (2..150) of |ratio - 1| is %f\n", sum);
 115
 116   TLine* line = new TLine(0, 1, 200, 1);
 117   line->SetLineWidth(2);
 118   line->Draw();
 119
 120   line = new TLine(0, 1.1, 200, 1.1);
 121   line->SetLineWidth(2);
 122   line->SetLineStyle(2);
 123   line->Draw();
 124   line = new TLine(0, 0.9, 200, 0.9);
 125   line->SetLineWidth(2);
 126   line->SetLineStyle(2);
 127   line->Draw();
 128
 129   canvas->Modified();
 130
 131   canvas->SaveAs(epsName);
 132
 133   return canvas;
 134 }
 135
 136 TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
 137 {
 138   // draws the 3 plots in the upper plot
 139   // draws the ratio between result1 and result2 in the lower plot
 140
 141   // normalize unfolded result to mc hist
 142   result1->Scale(1.0 / result1->Integral(2, 200));
 143   result1->Scale(mcHist->Integral(2, 200));
 144   result2->Scale(1.0 / result2->Integral(2, 200));
 145   result2->Scale(mcHist->Integral(2, 200));
 146
 147   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
 148   canvas->Range(0, 0, 1, 1);
 149
 150   TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
 151   pad1->Draw();
 152
 153   TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
 154   pad2->Draw();
 155
 156   pad1->SetRightMargin(0.05);
 157   pad2->SetRightMargin(0.05);
 158
 159   // no border between them
 160   pad1->SetBottomMargin(0);
 161   pad2->SetTopMargin(0);
 162
 163   pad1->cd();
 164
 165   mcHist->GetXaxis()->SetLabelSize(0.06);
 166   mcHist->GetYaxis()->SetLabelSize(0.06);
 167   mcHist->GetXaxis()->SetTitleSize(0.06);
 168   mcHist->GetYaxis()->SetTitleSize(0.06);
 169   mcHist->GetYaxis()->SetTitleOffset(0.6);
 170
 171   mcHist->GetXaxis()->SetRangeUser(0, 150);
 172
 173   mcHist->SetTitle(";true multiplicity;Entries");
 174   mcHist->SetStats(kFALSE);
 175
 176   mcHist->DrawCopy("HIST E");
 177   gPad->SetLogy();
 178
 179   result1->SetLineColor(2);
 180   result1->DrawCopy("SAME HISTE");
 181
 182   result2->SetLineColor(4);
 183   result2->DrawCopy("SAME HISTE");
 184
 185   TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
 186   legend->AddEntry(mcHist, "true distribution");
 187   legend->AddEntry(result1, "unfolded distribution (syst)");
 188   legend->AddEntry(result2, "unfolded distribution (normal)");
 189   legend->SetFillColor(0);
 190   legend->Draw();
 191
 192   pad2->cd();
 193   pad2->SetBottomMargin(0.15);
 194
 195   result1->GetXaxis()->SetLabelSize(0.06);
 196   result1->GetYaxis()->SetLabelSize(0.06);
 197   result1->GetXaxis()->SetTitleSize(0.06);
 198   result1->GetYaxis()->SetTitleSize(0.06);
 199   result1->GetYaxis()->SetTitleOffset(0.6);
 200
 201   result1->GetXaxis()->SetRangeUser(0, 150);
 202
 203   result1->SetTitle(";true multiplicity;Entries");
 204   result1->SetStats(kFALSE);
 205
 206   // calculate ratio
 207   result1->Sumw2();
 208   TH1* ratio = (TH1*) result1->Clone("ratio");
 209   result2->Sumw2();
 210   ratio->Divide(ratio, result2, 1, 1, "");
 211   ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
 212   ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
 213
 214   ratio->DrawCopy();
 215
 216   // get average of ratio
 217   Float_t sum = 0;
 218   for (Int_t i=2; i<=150; ++i)
 219   {
 220     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
 221   }
 222   sum /= 149;
 223
 224   printf("Average (2..150) of |ratio - 1| is %f\n", sum);
 225
 226   TLine* line = new TLine(0, 1, 150, 1);
 227   line->SetLineWidth(2);
 228   line->Draw();
 229
 230   line = new TLine(0, 1.1, 150, 1.1);
 231   line->SetLineWidth(2);
 232   line->SetLineStyle(2);
 233   line->Draw();
 234   line = new TLine(0, 0.9, 150, 0.9);
 235   line->SetLineWidth(2);
 236   line->SetLineStyle(2);
 237   line->Draw();
 238
 239   canvas->Modified();
 240
 241   canvas->SaveAs(epsName);
 242
 243   return canvas;
 244 }
 245
 246 TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE)
 247 {
 248   // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
 249   // a systematic effect
 250
 251   // normalize results
 252   result->Scale(1.0 / result->Integral(2, 200));
 253
 254   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
 255
 256   result->GetXaxis()->SetRangeUser(0, 150);
 257   result->SetStats(kFALSE);
 258
 259   for (Int_t n=0; n<nResultSyst; ++n)
 260   {
 261     resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
 262
 263     // calculate ratio
 264     TH1* ratio = (TH1*) result->Clone("ratio");
 265     ratio->Divide(ratio, resultSyst[n], 1, 1, "B");
 266     ratio->SetTitle(";true multiplicity;Ratio");
 267     ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
 268
 269     if (firstMarker)
 270       ratio->SetMarkerStyle(5);
 271
 272     ratio->SetLineColor(n+1);
 273
 274     ratio->DrawCopy((n == 0) ? ((firstMarker) ? "P" : "HIST") : "SAME HIST");
 275
 276     // get average of ratio
 277     Float_t sum = 0;
 278     for (Int_t i=2; i<=150; ++i)
 279       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
 280     sum /= 149;
 281
 282     printf("%d) Average (2..150) of |ratio - 1| is %f\n", n, sum);
 283   }
 284
 285   TLine* line = new TLine(0, 1, 150, 1);
 286   line->SetLineWidth(2);
 287   line->Draw();
 288
 289   line = new TLine(0, 1.1, 150, 1.1);
 290   line->SetLineWidth(2);
 291   line->SetLineStyle(2);
 292   line->Draw();
 293   line = new TLine(0, 0.9, 150, 0.9);
 294   line->SetLineWidth(2);
 295   line->SetLineStyle(2);
 296   line->Draw();
 297
 298   canvas->Modified();
 299
 300   canvas->SaveAs(epsName);
 301   canvas->SaveAs(Form("%s.gif", epsName.Data()));
 302
 303   return canvas;
 304 }
 305
 306 TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
 307 {
 308   // draws the ratios of each mc to the corresponding result
 309
 310   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
 311
 312   for (Int_t n=0; n<nResultSyst; ++n)
 313   {
 314     // normalize
 315     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
 316     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
 317
 318     result[n]->GetXaxis()->SetRangeUser(0, 150);
 319     result[n]->SetStats(kFALSE);
 320
 321     // calculate ratio
 322     TH1* ratio = (TH1*) result[n]->Clone("ratio");
 323     ratio->Divide(mc[n], ratio, 1, 1, "B");
 324     ratio->SetTitle(";true multiplicity;Ratio (true / unfolded)");
 325     ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
 326
 327     ratio->SetLineColor(n+1);
 328
 329     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
 330
 331     // get average of ratio
 332     Float_t sum = 0;
 333     for (Int_t i=2; i<=150; ++i)
 334       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
 335     sum /= 149;
 336
 337     printf("%d) Average (2..150) of |ratio - 1| is %f\n", n, sum);
 338   }
 339
 340   TLine* line = new TLine(0, 1, 150, 1);
 341   line->SetLineWidth(2);
 342   line->Draw();
 343
 344   line = new TLine(0, 1.1, 150, 1.1);
 345   line->SetLineWidth(2);
 346   line->SetLineStyle(2);
 347   line->Draw();
 348   line = new TLine(0, 0.9, 150, 0.9);
 349   line->SetLineWidth(2);
 350   line->SetLineStyle(2);
 351   line->Draw();
 352
 353   canvas->Modified();
 354
 355   canvas->SaveAs(epsName);
 356   canvas->SaveAs(Form("%s.gif", epsName.Data()));
 357
 358   return canvas;
 359 }
 360
 361 TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
 362 {
 363   // draws the ratios of each mc to the corresponding result
 364   // deducts from each ratio the ratio of mcBase / resultBase
 365
 366   // normalize
 367   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
 368   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
 369
 370   // calculate ratio
 371   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
 372   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
 373
 374   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
 375   canvas->SetRightMargin(0.05);
 376   canvas->SetTopMargin(0.05);
 377
 378   for (Int_t n=0; n<nResultSyst; ++n)
 379   {
 380     // normalize
 381     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
 382     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
 383
 384     result[n]->GetXaxis()->SetRangeUser(0, 150);
 385     result[n]->SetStats(kFALSE);
 386
 387     // calculate ratio
 388     TH1* ratio = (TH1*) result[n]->Clone("ratio");
 389     ratio->Divide(mc[n], ratio, 1, 1, "B");
 390     ratio->Add(ratioBase, -1);
 391
 392     ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
 393     ratio->GetYaxis()->SetRangeUser(-1, 1);
 394     ratio->SetLineColor(n+1);
 395     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
 396
 397     // get average of ratio
 398     Float_t sum = 0;
 399     for (Int_t i=2; i<=150; ++i)
 400       sum += TMath::Abs(ratio->GetBinContent(i));
 401     sum /= 149;
 402
 403     printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
 404   }
 405
 406   TLine* line = new TLine(0, 0, 150, 0);
 407   line->SetLineWidth(2);
 408   line->Draw();
 409
 410   line = new TLine(0, 0.1, 150, 0.1);
 411   line->SetLineWidth(2);
 412   line->SetLineStyle(2);
 413   line->Draw();
 414   line = new TLine(0, -0.1, 150, -0.1);
 415   line->SetLineWidth(2);
 416   line->SetLineStyle(2);
 417   line->Draw();
 418
 419   canvas->Modified();
 420
 421   canvas->SaveAs(epsName);
 422   canvas->SaveAs(Form("%s.gif", epsName.Data()));
 423
 424   return canvas;
 425 }
 426
 427 void Smooth(TH1* hist, Int_t windowWidth = 20)
 428 {
 429   TH1* clone = (TH1*) hist->Clone("clone");
 430   for (Int_t bin=3; bin<=clone->GetNbinsX(); ++bin)
 431   {
 432     Int_t min = TMath::Max(3, bin-windowWidth);
 433     Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
 434     Float_t average = clone->Integral(min, max) / (max - min + 1);
 435
 436     hist->SetBinContent(bin, average);
 437     hist->SetBinError(bin, 0);
 438   }
 439
 440   delete clone;
 441 }
 442
 443 TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
 444 {
 445   // draws the ratios of each mc to the corresponding result
 446   // deducts from each ratio the ratio of mcBase / resultBase
 447   // smoothens the ratios by a sliding window
 448
 449   // normalize
 450   resultBase->Scale(1.0 / resultBase->Integral(2, 200));
 451   mcBase->Scale(1.0 / mcBase->Integral(2, 200));
 452
 453   // calculate ratio
 454   TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
 455   ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
 456
 457   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
 458   canvas->SetRightMargin(0.05);
 459   canvas->SetTopMargin(0.05);
 460
 461   for (Int_t n=0; n<nResultSyst; ++n)
 462   {
 463     // normalize
 464     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
 465     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
 466
 467     result[n]->GetXaxis()->SetRangeUser(0, 150);
 468     result[n]->SetStats(kFALSE);
 469
 470     // calculate ratio
 471     TH1* ratio = (TH1*) result[n]->Clone("ratio");
 472     ratio->Divide(mc[n], ratio, 1, 1, "B");
 473     ratio->Add(ratioBase, -1);
 474
 475     new TCanvas; ratio->DrawCopy();
 476     Smooth(ratio);
 477     ratio->SetLineColor(1);
 478     ratio->DrawCopy("SAME");
 479
 480     canvas->cd();
 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<=150; ++i)
 489       sum += TMath::Abs(ratio->GetBinContent(i));
 490     sum /= 149;
 491
 492     printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
 493   }
 494
 495   TLine* line = new TLine(0, 0, 150, 0);
 496   line->SetLineWidth(2);
 497   line->Draw();
 498
 499   line = new TLine(0, 0.1, 150, 0.1);
 500   line->SetLineWidth(2);
 501   line->SetLineStyle(2);
 502   line->Draw();
 503   line = new TLine(0, -0.1, 150, -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 void DrawResiduals(TH1* measured, TH1* unfoldedFolded, const char* epsName)
 517 {
 518   // normalize
 519   unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, 200));
 520   unfoldedFolded->Scale(measured->Integral(2, 200));
 521
 522   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
 523   canvas->Range(0, 0, 1, 1);
 524
 525   TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
 526   pad1->Draw();
 527   pad1->SetGridx();
 528   pad1->SetGridy();
 529
 530   TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
 531   pad2->Draw();
 532   pad2->SetGridx();
 533   pad2->SetGridy();
 534
 535   TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
 536   pad3->SetGridx();
 537   pad3->SetGridy();
 538   pad3->SetRightMargin(0.05);
 539   pad3->SetTopMargin(0.05);
 540   pad3->Draw();
 541
 542   pad1->SetRightMargin(0.05);
 543   pad2->SetRightMargin(0.05);
 544
 545   // no border between them
 546   pad1->SetBottomMargin(0);
 547   pad2->SetTopMargin(0);
 548
 549   pad1->cd();
 550
 551   measured->GetXaxis()->SetLabelSize(0.06);
 552   measured->GetYaxis()->SetLabelSize(0.06);
 553   measured->GetXaxis()->SetTitleSize(0.06);
 554   measured->GetYaxis()->SetTitleSize(0.06);
 555   measured->GetYaxis()->SetTitleOffset(0.6);
 556
 557   measured->GetXaxis()->SetRangeUser(0, 150);
 558
 559   measured->SetTitle(";measured multiplicity;Entries");
 560   measured->SetStats(kFALSE);
 561
 562   measured->DrawCopy("HIST");
 563   gPad->SetLogy();
 564
 565   unfoldedFolded->SetMarkerStyle(5);
 566   unfoldedFolded->SetMarkerColor(2);
 567   unfoldedFolded->SetLineColor(0);
 568   unfoldedFolded->DrawCopy("SAME P");
 569
 570   TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
 571   legend->AddEntry(measured, "measured distribution");
 572   legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution");
 573   legend->SetFillColor(0);
 574   legend->Draw();
 575
 576   pad2->cd();
 577   pad2->SetBottomMargin(0.15);
 578
 579   // calculate ratio
 580   measured->Sumw2();
 581   TH1* residual = (TH1*) measured->Clone("residual");
 582   unfoldedFolded->Sumw2();
 583
 584   residual->Add(unfoldedFolded, -1);
 585
 586   // projection
 587   TH1* residualHist = new TH1F("residualHist", ";", 15, -3, 3);
 588
 589   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
 590   {
 591     if (measured->GetBinError(i) > 0)
 592     {
 593       residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
 594       residual->SetBinError(i, 1);
 595
 596       residualHist->Fill(residual->GetBinContent(i));
 597     }
 598     else
 599     {
 600       residual->SetBinContent(i, 0);
 601       residual->SetBinError(i, 0);
 602     }
 603   }
 604
 605   residual->GetYaxis()->SetTitle("Residuals   1/e (M - R #otimes U)");
 606   residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
 607   residual->DrawCopy();
 608
 609   TLine* line = new TLine(-0.5, 0, 150.5, 0);
 610   line->SetLineWidth(2);
 611   line->Draw();
 612
 613   pad3->cd();
 614   residualHist->SetStats(kFALSE);
 615   residualHist->GetXaxis()->SetLabelSize(0.08);
 616   residualHist->Fit("gaus");
 617   residualHist->Draw();
 618
 619   canvas->Modified();
 620   canvas->SaveAs(canvas->GetName());
 621
 622   //const char* epsName2 = "proj.eps";
 623   //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
 624   //canvas->SetGridx();
 625   //canvas->SetGridy();
 626
 627   //canvas->SaveAs(canvas->GetName());
 628 }
 629
 630 void bayesianExample()
 631 {
 632   TStopwatch watch;
 633   watch.Start();
 634
 635   gSystem->Load("libPWG0base");
 636
 637   TFile::Open(correctionFile);
 638   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 639   mult->LoadHistograms("Multiplicity");
 640
 641   TFile::Open(measuredFile);
 642   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 643   mult2->LoadHistograms("Multiplicity");
 644
 645   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 646
 647   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
 648
 649   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 650   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 651
 652   //mult->DrawComparison("bayesianExample", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
 653   DrawResultRatio(mcHist, result, "bayesianExample.eps");
 654
 655   // draw residual plot
 656
 657   // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
 658   TH2* convoluted = mult->CalculateMultiplicityESD(result, etaRange);
 659   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
 660
 661   TH1* measured = mult2->GetMultiplicityESD(etaRange)->ProjectionY("measured");
 662
 663   DrawResiduals(measured, convolutedProj, "bayesianResiduals.eps");
 664
 665   watch.Stop();
 666   watch.Print();
 667 }
 668
 669 void chi2FluctuationResult()
 670 {
 671   gSystem->Load("libPWG0base");
 672
 673   TFile::Open(correctionFile);
 674   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 675   mult->LoadHistograms("Multiplicity");
 676
 677   TFile::Open(measuredFile);
 678   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 679   mult2->LoadHistograms("Multiplicity");
 680
 681   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 682   mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
 683   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 684
 685   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 686   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 687
 688   mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
 689
 690   TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_3");
 691   canvas->SaveAs("chi2FluctuationResult.eps");
 692 }
 693
 694 void chi2Example()
 695 {
 696   TStopwatch watch;
 697   watch.Start();
 698
 699   gSystem->Load("libPWG0base");
 700
 701   TFile::Open(correctionFile);
 702   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 703   mult->LoadHistograms("Multiplicity");
 704
 705   TFile::Open(measuredFile);
 706   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 707   mult2->LoadHistograms("Multiplicity");
 708
 709   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 710   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
 711   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 712
 713   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 714   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 715
 716   DrawResultRatio(mcHist, result, "chi2Example.eps");
 717
 718   watch.Stop();
 719   watch.Print();
 720 }
 721
 722 void chi2ExampleTPC()
 723 {
 724   TStopwatch watch;
 725   watch.Start();
 726
 727   gSystem->Load("libPWG0base");
 728
 729   TFile::Open(correctionFileTPC);
 730   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 731   mult->LoadHistograms("Multiplicity");
 732
 733   TFile::Open(measuredFileTPC);
 734   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 735   mult2->LoadHistograms("Multiplicity");
 736
 737   mult->SetMultiplicityESD(etaRangeTPC, mult2->GetMultiplicityESD(etaRangeTPC));
 738   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
 739   mult->ApplyMinuitFit(etaRangeTPC, kFALSE, AliMultiplicityCorrection::kTrVtx);
 740
 741   TH1* mcHist = mult2->GetMultiplicityVtx(etaRangeTPC)->ProjectionY("mymc");
 742   TH1* result = mult->GetMultiplicityESDCorrected(etaRangeTPC);
 743
 744   DrawResultRatio(mcHist, result, "chi2ExampleTPC.eps");
 745
 746   watch.Stop();
 747   watch.Print();
 748 }
 749
 750 void bayesianNBD()
 751 {
 752   gSystem->Load("libPWG0base");
 753   TFile::Open("multiplicityMC_3M.root");
 754
 755   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 756   mult->LoadHistograms("Multiplicity");
 757
 758   TFile::Open("multiplicityMC_3M_NBD.root");
 759   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 760   mult2->LoadHistograms("Multiplicity");
 761
 762   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 763   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1.0, 100);
 764
 765   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 766
 767   mcHist->Sumw2();
 768   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 769
 770   //mult->DrawComparison("bayesianNBD", etaRange, kFALSE, kTRUE, mcHist);
 771   DrawResultRatio(mcHist, result, "bayesianNBD.eps");
 772 }
 773
 774 void minimizationNBD()
 775 {
 776   gSystem->Load("libPWG0base");
 777   TFile::Open("multiplicityMC_3M.root");
 778
 779   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 780   mult->LoadHistograms("Multiplicity");
 781
 782   TFile::Open("multiplicityMC_3M_NBD.root");
 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
 792   mcHist->Sumw2();
 793   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 794
 795   //mult->DrawComparison("minimizationNBD", etaRange, kFALSE, kTRUE, mcHist);
 796   DrawResultRatio(mcHist, result, "minimizationNBD.eps");
 797 }
 798
 799 void minimizationInfluenceAlpha()
 800 {
 801   gSystem->Load("libPWG0base");
 802
 803   TFile::Open(measuredFile);
 804   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
 805   mult2->LoadHistograms("Multiplicity");
 806
 807   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 808   mcHist->Scale(1.0 / mcHist->Integral());
 809   mcHist->GetXaxis()->SetRangeUser(0, 200);
 810   mcHist->SetStats(kFALSE);
 811   mcHist->SetTitle(";true multiplicity;P_{N}");
 812
 813   TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 300);
 814   canvas->Divide(3, 1);
 815
 816   TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
 817
 818   TH1* hist1 = gFile->Get("MinuitChi2_00_2_100.000000");
 819   TH1* hist2 = gFile->Get("MinuitChi2_03_2_100000.000000");
 820   TH1* hist3 = gFile->Get("MinuitChi2_06_2_100000000.000000");
 821
 822   mcHist->Rebin(2);  mcHist->Scale(0.5);
 823   hist1->Rebin(2);   hist1->Scale(0.5);
 824   hist2->Rebin(2);   hist2->Scale(0.5);
 825   hist3->Rebin(2);   hist3->Scale(0.5);
 826
 827   mcHist->GetXaxis()->SetRangeUser(0, 200);
 828
 829   canvas->cd(1);
 830   gPad->SetLogy();
 831   mcHist->Draw();
 832   hist1->SetMarkerStyle(5);
 833   hist1->SetMarkerColor(2);
 834   hist1->Draw("SAME PE");
 835
 836   canvas->cd(2);
 837   gPad->SetLogy();
 838   mcHist->Draw();
 839   hist2->SetMarkerStyle(5);
 840   hist2->SetMarkerColor(2);
 841   hist2->Draw("SAME PE");
 842
 843   canvas->cd(3);
 844   gPad->SetLogy();
 845   mcHist->Draw();
 846   hist3->SetMarkerStyle(5);
 847   hist3->SetMarkerColor(2);
 848   hist3->Draw("SAME PE");
 849
 850   canvas->SaveAs("minimizationInfluenceAlpha.eps");
 851 }
 852
 853 void NBDFit()
 854 {
 855   gSystem->Load("libPWG0base");
 856
 857   TFile::Open(correctionFile);
 858   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 859   mult->LoadHistograms("Multiplicity");
 860
 861   TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
 862   fCurrentESD->Sumw2();
 863   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 864
 865   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])");
 866   func->SetParNames("scaling", "averagen", "k");
 867   func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
 868   func->SetParLimits(1, 0.001, 1000);
 869   func->SetParLimits(2, 0.001, 1000);
 870   func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
 871
 872   TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
 873   lognormal->SetParNames("scaling", "mean", "sigma");
 874   lognormal->SetParameters(1, 1, 1);
 875   lognormal->SetParLimits(0, 0, 10);
 876   lognormal->SetParLimits(1, 0, 100);
 877   lognormal->SetParLimits(2, 1e-3, 10);
 878
 879   TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
 880   fCurrentESD->SetStats(kFALSE);
 881   fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
 882   fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
 883   fCurrentESD->Draw("HIST");
 884   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
 885   fCurrentESD->Fit(func, "W0", "", 0, 50);
 886   func->SetRange(0, 100);
 887   func->Draw("SAME");
 888   printf("chi2 = %f\n", func->GetChisquare());
 889
 890   fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
 891   lognormal->SetLineColor(2);
 892   lognormal->SetLineStyle(2);
 893   lognormal->SetRange(0, 100);
 894   lognormal->Draw("SAME");
 895
 896   canvas->SaveAs("NBDFit.eps");
 897 }
 898
 899 void DifferentSamples()
 900 {
 901   // data generated by runMultiplicitySelector.C DifferentSamples
 902
 903   const char* name = "DifferentSamples";
 904
 905   TFile* file = TFile::Open(Form("%s.root", name));
 906
 907   TCanvas* canvas = new TCanvas(name, name, 800, 600);
 908   canvas->Divide(2, 2);
 909
 910   for (Int_t i=0; i<4; ++i)
 911   {
 912     canvas->cd(i+1);
 913     gPad->SetTopMargin(0.05);
 914     gPad->SetRightMargin(0.05);
 915     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", i));
 916     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", i));
 917     TH1* mc = (TH1*) file->Get(Form("mc_%d", i));
 918     mc->Sumw2();
 919
 920     chi2Result->Divide(chi2Result, mc, 1, 1, "");
 921     bayesResult->Divide(bayesResult, mc, 1, 1, "");
 922
 923     chi2Result->SetTitle(";true multiplicity;unfolded measured/MC");
 924     chi2Result->GetXaxis()->SetRangeUser(0, 150);
 925     chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
 926     chi2Result->GetYaxis()->SetTitleOffset(1.2);
 927     chi2Result->SetLineColor(1);
 928     chi2Result->SetStats(kFALSE);
 929
 930     bayesResult->SetStats(kFALSE);
 931     bayesResult->SetLineColor(2);
 932
 933     chi2Result->DrawCopy("HIST");
 934     bayesResult->DrawCopy("SAME HIST");
 935
 936     TLine* line = new TLine(0, 1, 150, 1);
 937     line->Draw();
 938   }
 939
 940   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 941 }
 942
 943 void StartingConditions()
 944 {
 945   // data generated by runMultiplicitySelector.C StartingConditions
 946
 947   const char* name = "StartingConditions";
 948
 949   TFile* file = TFile::Open(Form("%s.root", name));
 950
 951   TCanvas* canvas = new TCanvas(name, name, 800, 400);
 952   canvas->Divide(2, 1);
 953
 954   TH1* mc = (TH1*) file->Get("mc");
 955   mc->Sumw2();
 956   mc->Scale(1.0 / mc->Integral());
 957
 958   Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
 959
 960   for (Int_t i=0; i<6; ++i)
 961   {
 962     Int_t id = i;
 963     if (id > 2)
 964       id += 2;
 965
 966     TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
 967     TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
 968
 969     chi2Result->Divide(chi2Result, mc, 1, 1, "");
 970     bayesResult->Divide(bayesResult, mc, 1, 1, "");
 971
 972     chi2Result->SetTitle("a) #chi^{2} minimization;true multiplicity;unfolded / MC");
 973     chi2Result->GetXaxis()->SetRangeUser(0, 150);
 974     chi2Result->GetYaxis()->SetRangeUser(0.8, 1.2);
 975     chi2Result->GetYaxis()->SetTitleOffset(1.5);
 976     //chi2Result->SetMarkerStyle(marker[i]);
 977     chi2Result->SetLineColor(i+1);
 978     chi2Result->SetStats(kFALSE);
 979
 980     bayesResult->SetTitle("b) Bayesian method;true multiplicity;unfolded / MC");
 981     bayesResult->GetXaxis()->SetRangeUser(0, 150);
 982     bayesResult->GetYaxis()->SetRangeUser(0.8, 1.2);
 983     bayesResult->GetYaxis()->SetTitleOffset(1.5);
 984     bayesResult->SetStats(kFALSE);
 985     bayesResult->SetLineColor(2);
 986     bayesResult->SetLineColor(i+1);
 987
 988     canvas->cd(1);
 989     gPad->SetLeftMargin(0.12);
 990     chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
 991
 992     canvas->cd(2);
 993     gPad->SetLeftMargin(0.12);
 994     bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
 995
 996     //TLine* line = new TLine(0, 1, 150, 1);
 997     //line->Draw();
 998   }
 999
1000   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1001 }
1002
1003 void StatisticsPlot()
1004 {
1005   const char* name = "StatisticsPlot";
1006
1007   TFile* file = TFile::Open(Form("%s.root", name));
1008
1009   TCanvas* canvas = new TCanvas(name, name, 600, 400);
1010
1011   TGraph* fitResultsChi2 = file->Get("fitResultsChi2");
1012   fitResultsChi2->SetTitle(";number of measured events;P_{1}");
1013   fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
1014   fitResultsChi2->Draw("AP");
1015
1016   TF1* f = new TF1("f", "[0]/x", 1, 1e4);
1017   fitResultsChi2->Fit(f);
1018
1019   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1020
1021   TH1* mc[5];
1022   TH1* result[5];
1023
1024   const char* plotname = "chi2Result";
1025
1026   const char* name = "StatisticsPlotRatios";
1027
1028   TCanvas* canvas = new TCanvas(name, name, 600, 400);
1029
1030   for (Int_t i=0; i<5; ++i)
1031   {
1032     mc[i] = (TH1*) file->Get(Form("mc_%d", i));
1033     result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
1034
1035     result[i]->SetLineColor(i+1);
1036     result[i]->Draw(((i == 0) ? "" : "SAME"));
1037   }
1038 }
1039
1040 void SystematicLowEfficiency()
1041 {
1042   gSystem->Load("libPWG0base");
1043
1044   TFile::Open(correctionFile);
1045   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1046   mult->LoadHistograms("Multiplicity");
1047
1048   // calculate result with systematic effect
1049   TFile::Open("multiplicityMC_100k_1_loweff.root");
1050   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1051   mult2->LoadHistograms("Multiplicity");
1052
1053   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1054   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1055   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1056
1057   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1058   TH1* result1 = mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
1059
1060   DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
1061
1062   // calculate normal result
1063   TFile::Open("multiplicityMC_100k_1.root");
1064   mult2->LoadHistograms("Multiplicity");
1065
1066   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1067   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1068
1069   TH1* result2 = mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1070
1071   DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
1072
1073   Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
1074 }
1075
1076 void SystematicMisalignment()
1077 {
1078   gSystem->Load("libPWG0base");
1079
1080   TFile::Open(correctionFile);
1081   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1082   mult->LoadHistograms("Multiplicity");
1083
1084   TFile::Open("multiplicityMC_100k_fullmis.root");
1085   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1086   mult2->LoadHistograms("Multiplicity");
1087
1088   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1089   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1090   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1091
1092   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1093   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1094
1095   DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
1096 }
1097
1098 void EfficiencySpecies(const char* fileName = "multiplicityMC_400k_syst.root")
1099 {
1100   gSystem->Load("libPWG0base");
1101
1102   TFile::Open(fileName);
1103   AliCorrection* correction[4];
1104
1105   TCanvas* canvas = new TCanvas("EfficiencySpeciesSPD", "EfficiencySpeciesSPD", 600, 400);
1106   canvas->SetGridx();
1107   canvas->SetGridy();
1108   canvas->SetRightMargin(0.05);
1109   canvas->SetTopMargin(0.05);
1110
1111   TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
1112   legend->SetFillColor(0);
1113
1114   Int_t marker[] = {24, 25, 26};
1115   Int_t color[] = {1, 2, 4};
1116
1117   Float_t below = 0;
1118   Float_t total = 0;
1119
1120   for (Int_t i=0; i<3; ++i)
1121   {
1122     Printf("correction %d", i);
1123
1124     TString name; name.Form("correction_%d", i);
1125     correction[i] = new AliCorrection(name, name);
1126     correction[i]->LoadHistograms();
1127
1128     TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
1129     TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
1130
1131     // limit vtx axis
1132     gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
1133     meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
1134
1135     // 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)
1136     /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1137       for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
1138       {
1139         gene->SetBinContent(x, 0, z, 0);
1140         gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1141         meas->SetBinContent(x, 0, z, 0);
1142         meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
1143       }*/
1144
1145     // limit eta axis
1146     gene->GetYaxis()->SetRangeUser(-0.49, 0.49);
1147     meas->GetYaxis()->SetRangeUser(-0.49, 0.49);
1148
1149     TH1* genePt = gene->Project3D(Form("z_%d", i));
1150     TH1* measPt = meas->Project3D(Form("z_%d", i));
1151
1152     genePt->Sumw2();
1153     measPt->Sumw2();
1154
1155     effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
1156     effPt->Reset();
1157     effPt->Divide(measPt, genePt, 1, 1, "B");
1158
1159     Int_t bin = 0;
1160     for (bin=20; bin>=1; bin--)
1161     {
1162       if (effPt->GetBinContent(bin) < 0.5)
1163         break;
1164     }
1165
1166     Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
1167
1168     Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
1169     Printf("%.4f of the particles are below that momentum", fraction);
1170
1171     below += genePt->Integral(1, bin);
1172     total += genePt->Integral();
1173
1174     effPt->SetLineColor(color[i]);
1175     effPt->SetMarkerColor(color[i]);
1176     effPt->SetMarkerStyle(marker[i]);
1177
1178     effPt->GetXaxis()->SetRangeUser(0, 1);
1179     effPt->GetYaxis()->SetRangeUser(0, 1);
1180
1181     effPt->SetStats(kFALSE);
1182     effPt->SetTitle("");
1183     effPt->GetYaxis()->SetTitle("Efficiency");
1184
1185     effPt->DrawCopy((i == 0) ? "" : "SAME");
1186
1187     legend->AddEntry(effPt, ((i == 0) ? "#pi" : ((i == 1) ? "K" : "p")));
1188   }
1189
1190   Printf("In total %.4f of the particles are below their pt cut off", (Float_t) below / total);
1191
1192   legend->Draw();
1193
1194   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1195 }
1196
1197 void ParticleSpeciesComparison1()
1198 {
1199   gSystem->Load("libPWG0base");
1200
1201   const char* fileNameMC = "multiplicityMC_400k_syst_species.root";
1202   const char* fileNameESD = "multiplicityMC_100k_syst.root";
1203   Bool_t chi2 = 0;
1204
1205   TFile::Open(fileNameESD);
1206   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
1207   TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
1208
1209   TH1* results[10];
1210
1211   // loop over cases (normal, enhanced/reduced ratios)
1212   Int_t nMax = 7;
1213   for (Int_t i = 0; i<nMax; ++i)
1214   {
1215     TString folder;
1216     folder.Form("Multiplicity_%d", i);
1217
1218     AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
1219
1220     TFile::Open(fileNameMC);
1221     mult->LoadHistograms();
1222
1223     mult->SetMultiplicityESD(etaRange, hist);
1224
1225     if (chi2)
1226     {
1227       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1228       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1229       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1230     }
1231     else
1232     {
1233       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1234       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1235     }
1236
1237     //Float_t averageRatio = 0;
1238     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1239
1240     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1241
1242     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1243   }
1244
1245   DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
1246
1247   TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
1248
1249   for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
1250   {
1251     results[0]->SetBinError(i, 0);
1252     mc->SetBinError(i, 0);
1253   }
1254
1255   TCanvas* canvas = DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps");
1256
1257   TFile::Open("bayesianUncertainty_400k_100k_syst.root");
1258   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1259   errorHist->SetLineColor(1);
1260   errorHist->SetLineWidth(2);
1261   TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1262   for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1263   {
1264     errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1265     errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1266   }
1267
1268   errorHist->DrawCopy("SAME");
1269   errorHist2->DrawCopy("SAME");
1270
1271   canvas->SaveAs(canvas->GetName());
1272
1273   TCanvas* canvas2 = DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE);
1274
1275   errorHist->DrawCopy("SAME");
1276   errorHist2->DrawCopy("SAME");
1277
1278   canvas2->SaveAs(canvas2->GetName());
1279 }
1280
1281 void ParticleSpeciesComparison2()
1282 {
1283   gSystem->Load("libPWG0base");
1284
1285   const char* fileNameMC = "multiplicityMC_400k_syst.root";
1286   const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
1287   Bool_t chi2 = 0;
1288
1289   TFile::Open(fileNameMC);
1290   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1291   mult->LoadHistograms();
1292
1293   TH1* mc[10];
1294   TH1* results[10];
1295
1296   // loop over cases (normal, enhanced/reduced ratios)
1297   Int_t nMax = 7;
1298   for (Int_t i = 0; i<nMax; ++i)
1299   {
1300     TString folder;
1301     folder.Form("Multiplicity_%d", i);
1302
1303     AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
1304
1305     TFile::Open(fileNameESD);
1306     mult2->LoadHistograms();
1307
1308     mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1309
1310     if (chi2)
1311     {
1312       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
1313       mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
1314       //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
1315     }
1316     else
1317     {
1318       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1319       //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
1320     }
1321
1322     //Float_t averageRatio = 0;
1323     //mult->GetComparisonResults(0, 0, 0, &averageRatio);
1324
1325     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
1326
1327     TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
1328     mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
1329
1330     //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
1331     //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
1332
1333     //Printf("Case %d. Average ratio is %f", i, averageRatio);
1334   }
1335
1336   DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
1337 }
1338
1339 TH1* Invert(TH1* eff)
1340 {
1341   // calculate corr = 1 / eff
1342
1343   TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
1344   corr->Reset();
1345
1346   for (Int_t i=1; i<=eff->GetNbinsX(); i++)
1347   {
1348     if (eff->GetBinContent(i) > 0)
1349     {
1350       corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
1351       corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
1352     }
1353   }
1354
1355   return corr;
1356 }
1357
1358 void TriggerVertexCorrection()
1359 {
1360   //
1361   // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
1362   //
1363
1364   gSystem->Load("libPWG0base");
1365
1366   TFile::Open(correctionFile);
1367   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1368   mult->LoadHistograms("Multiplicity");
1369
1370   TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
1371   TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
1372
1373   TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
1374
1375   corrINEL->SetStats(kFALSE);
1376   corrINEL->GetXaxis()->SetRangeUser(0, 30);
1377   corrINEL->GetYaxis()->SetRangeUser(0, 5);
1378   corrINEL->SetTitle(";true multiplicity;correction factor");
1379   corrINEL->SetMarkerStyle(22);
1380   corrINEL->Draw("PE");
1381
1382   corrMB->SetStats(kFALSE);
1383   corrMB->SetLineColor(2);
1384   corrMB->SetMarkerStyle(25);
1385   corrMB->SetMarkerColor(2);
1386   corrMB->Draw("SAME PE");
1387
1388   TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
1389   legend->SetFillColor(0);
1390   legend->AddEntry(corrINEL, "correction to inelastic sample");
1391   legend->AddEntry(corrMB, "correction to minimum bias sample");
1392
1393   legend->Draw();
1394
1395   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1396 }
1397
1398 void bayesianUncertainty(Bool_t mc = kFALSE)
1399 {
1400   gSystem->Load("libPWG0base");
1401
1402   TFile::Open(correctionFile);
1403   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1404   mult->LoadHistograms("Multiplicity");
1405
1406   TFile::Open(measuredFile);
1407   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1408   mult2->LoadHistograms("Multiplicity");
1409
1410   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1411
1412   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1413
1414   TH1* errorResponse = mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, 1, 100, ((mc) ? mcHist : 0))->Clone("errorResponse");
1415   TH1* errorMeasured = mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, 1, 100, ((mc) ? mcHist : 0))->Clone("errorMeasured");
1416   TH1* errorBoth = mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0))->Clone("errorBoth");
1417
1418   if (!mc)
1419   {
1420     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1421     DrawResultRatio(mcHist, result, "bayesianUncertainty2.eps");
1422   }
1423
1424   TCanvas* canvas = new TCanvas("bayesianUncertainty", "bayesianUncertainty", 600, 400);
1425   canvas->SetGridx();
1426   canvas->SetGridy();
1427   canvas->SetRightMargin(0.05);
1428   canvas->SetTopMargin(0.05);
1429
1430   errorResponse->SetLineColor(1);
1431   errorResponse->GetXaxis()->SetRangeUser(0, 200);
1432   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
1433   errorResponse->SetStats(kFALSE);
1434   errorResponse->SetTitle(";true multiplicity;Uncertainty");
1435
1436   errorResponse->Draw();
1437
1438   errorMeasured->SetLineColor(2);
1439   errorMeasured->Draw("SAME");
1440
1441   errorBoth->SetLineColor(3);
1442   errorBoth->Draw("SAME");
1443
1444   Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
1445   Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
1446   Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
1447
1448   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1449
1450   TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
1451   errorResponse->Write();
1452   errorMeasured->Write();
1453   errorBoth->Write();
1454   file->Close();
1455 }
1456
1457 void EfficiencyComparison(Int_t eventType = 2)
1458 {
1459  const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
1460
1461   gSystem->Load("libPWG0base");
1462
1463   TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
1464   canvas->SetGridx();
1465   canvas->SetGridy();
1466   canvas->SetRightMargin(0.05);
1467   canvas->SetTopMargin(0.05);
1468
1469   AliMultiplicityCorrection* data[4];
1470   TH1* effArray[4];
1471
1472   Int_t markers[] = { 24, 25, 26, 5 };
1473   Int_t colors[] = { 1, 2, 3, 4 };
1474
1475   TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
1476   legend->SetFillColor(0);
1477
1478   TH1* effError = 0;
1479
1480   for (Int_t i=0; i<4; ++i)
1481   {
1482     TString name;
1483     name.Form("Multiplicity_%d", i);
1484
1485     TFile::Open(files[i]);
1486     data[i] = new AliMultiplicityCorrection(name, name);
1487
1488     if (i < 3)
1489     {
1490       data[i]->LoadHistograms("Multiplicity");
1491     }
1492     else
1493       data[i]->LoadHistograms("Multiplicity_0");
1494
1495     TH1* eff = (TH1*) data[i]->GetEfficiency(3, eventType)->Clone(Form("eff_%d", i));
1496     effArray[i] = eff;
1497
1498     eff->GetXaxis()->SetRangeUser(0, 15);
1499     eff->GetYaxis()->SetRangeUser(0, 1.1);
1500     eff->SetStats(kFALSE);
1501     eff->SetTitle(";true multiplicity;Efficiency");
1502     eff->SetLineColor(colors[i]);
1503     eff->SetMarkerColor(colors[i]);
1504     eff->SetMarkerStyle(markers[i]);
1505
1506     if (i == 3)
1507     {
1508       for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1509         eff->SetBinError(bin, 0);
1510
1511       // loop over cross section combinations
1512       for (Int_t j=1; j<7; ++j)
1513       {
1514         AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
1515         mult->LoadHistograms(Form("Multiplicity_%d", j));
1516
1517         TH1* eff2 = mult->GetEfficiency(3, eventType);
1518
1519         for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
1520         {
1521           // TODO we could also do assymetric errors here
1522           Float_t deviation = 10.0 * TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
1523
1524           eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), deviation));
1525         }
1526       }
1527
1528       for (Int_t bin=1; bin<=20; bin++)
1529         if (eff->GetBinContent(bin) > 0)
1530           Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
1531
1532       effError = (TH1*) eff2->Clone("effError");
1533       effError->Reset();
1534
1535       for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
1536         if (eff->GetBinContent(bin) > 0)
1537           effError->SetBinContent(bin, eff->GetBinError(bin) / eff->GetBinContent(bin));
1538
1539       effError->DrawCopy("SAME HIST");
1540     }
1541
1542     eff->SetBinContent(1, 0);
1543     eff->SetBinError(1, 0);
1544
1545     canvas->cd();
1546     if (i == 0)
1547     {
1548       eff->DrawCopy("P");
1549     }
1550     else
1551       eff->DrawCopy("SAME P");
1552
1553     legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
1554   }
1555
1556   legend->AddEntry(effError, "relative syst. uncertainty #times 10");
1557
1558   legend->Draw();
1559
1560   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1561 }
1562
1563 void ModelDependencyPlot()
1564 {
1565   gSystem->Load("libPWG0base");
1566
1567   TFile::Open("multiplicityMC_3M.root");
1568   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1569   mult->LoadHistograms("Multiplicity");
1570
1571   TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
1572
1573   TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
1574   canvas->SetGridx();
1575   canvas->SetGridy();
1576   //canvas->SetRightMargin(0.05);
1577   //canvas->SetTopMargin(0.05);
1578
1579   canvas->Divide(2, 1);
1580
1581   canvas->cd(2);
1582   gPad->SetLogy();
1583
1584   Int_t selectedMult = 30;
1585   Int_t yMax = 200000;
1586
1587   TH1* full = proj->ProjectionX("full");
1588   TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult));
1589
1590   full->SetStats(kFALSE);
1591   full->GetXaxis()->SetRangeUser(0, 200);
1592   full->GetYaxis()->SetRangeUser(5, yMax);
1593   full->SetTitle(";multiplicity");
1594
1595   selected->SetLineColor(0);
1596   selected->SetMarkerColor(2);
1597   selected->SetMarkerStyle(7);
1598
1599   full->Draw();
1600   selected->Draw("SAME P");
1601
1602   TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
1603   legend->SetFillColor(0);
1604   legend->AddEntry(full, "true");
1605   legend->AddEntry(selected, "measured");
1606   legend->Draw();
1607
1608   TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
1609   line->SetLineWidth(2);
1610   line->Draw();
1611
1612   canvas->cd(1);
1613   gPad->SetLogy();
1614
1615   TH1* full = proj->ProjectionY("full2");
1616   TH1* selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
1617
1618   full->SetStats(kFALSE);
1619   full->GetXaxis()->SetRangeUser(0, 200);
1620   full->GetYaxis()->SetRangeUser(5, yMax);
1621   full->SetTitle(";multiplicity");
1622
1623   full->SetLineColor(0);
1624   full->SetMarkerColor(2);
1625   full->SetMarkerStyle(7);
1626
1627   full->Draw("P");
1628   selected->Draw("SAME");
1629
1630   legend->Draw();
1631
1632   TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
1633   line->SetLineWidth(2);
1634   line->Draw();
1635
1636   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1637 }
1638
1639 void SystematicpTSpectrum()
1640 {
1641   gSystem->Load("libPWG0base");
1642
1643   TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
1644   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1645   mult->LoadHistograms("Multiplicity");
1646
1647   TFile::Open("multiplicityMC_100k_syst.root");
1648   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1649   mult2->LoadHistograms("Multiplicity");
1650
1651   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1652   mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1653   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1654
1655   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1656   TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
1657
1658   DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
1659 }
1660
1661 // to be deleted
1662 /*void covMatrix(Bool_t mc = kTRUE)
1663 {
1664   gSystem->Load("libPWG0base");
1665
1666   TFile::Open(correctionFile);
1667   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1668   mult->LoadHistograms("Multiplicity");
1669
1670   TFile::Open(measuredFile);
1671   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1672   mult2->LoadHistograms("Multiplicity");
1673
1674   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1675
1676   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
1677
1678   mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
1679 }*/
1680
1681 Double_t FitPtFunc(Double_t *x, Double_t *par)
1682 {
1683   Double_t xx = x[0];
1684
1685   Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
1686   Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
1687
1688   const Float_t kTransitionWidth = 0;
1689
1690   // power law part
1691   if (xx < par[0] - kTransitionWidth)
1692   {
1693     return val1;
1694   }
1695   /*else if (xx < par[0] + kTransitionWidth)
1696   {
1697     // smooth transition
1698     Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
1699     return (1 - factor) * val1 + factor * val2;
1700   }*/
1701   else
1702   {
1703     return val2;
1704   }
1705 }
1706
1707 void FitPt(const char* fileName = "firstplots100k_truept.root")
1708 {
1709   gSystem->Load("libPWG0base");
1710
1711   TFile::Open(fileName);
1712
1713   /*
1714   // merge corrections
1715   AliCorrection* correction[4];
1716   TList list;
1717
1718   for (Int_t i=0; i<4; ++i)
1719   {
1720     Printf("correction %d", i);
1721
1722     TString name; name.Form("correction_%d", i);
1723     correction[i] = new AliCorrection(name, name);
1724     correction[i]->LoadHistograms();
1725
1726     if (i > 0)
1727       list.Add(correction[i]);
1728   }
1729
1730   correction[0]->Merge(&list);
1731
1732   TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
1733
1734   // limit vtx, eta axis
1735   gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
1736   gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
1737
1738   TH1* genePt = gene->Project3D("z");*/
1739   TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
1740   genePt->Sumw2();
1741
1742   //genePt->Scale(1.0 / genePt->Integral());
1743
1744   // normalize by bin width
1745   for (Int_t x=1; x<genePt->GetNbinsX(); x++)
1746     genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
1747
1748   /// genePt->GetXaxis()->GetBinCenter(x));
1749
1750   genePt->GetXaxis()->SetRangeUser(0, 7.9);
1751   genePt->GetYaxis()->SetTitle("a.u.");
1752
1753   TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
1754   //func->SetLineColor(2);
1755   func->SetParameters(1, -1, 1, 1, 1);
1756   func->SetParLimits(3, 1, 10);
1757   func->SetParLimits(4, 0, 10);
1758
1759   //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
1760
1761   //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
1762   //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
1763   //func->FixParameter(0, 0.314);
1764   //func->SetParLimits(0, 0.1, 0.3);
1765
1766   genePt->SetMarkerStyle(25);
1767   genePt->SetTitle("");
1768   genePt->SetStats(kFALSE);
1769   genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
1770   //func->Draw("SAME");
1771
1772   // fit only exp. part
1773   func->SetParameters(1, -1);
1774   func->FixParameter(2, 0);
1775   func->FixParameter(3, 1);
1776   func->FixParameter(4, 1);
1777   genePt->Fit(func, "0", "", 0.2, 1);
1778
1779   new TCanvas;
1780   genePt->DrawCopy("P");
1781   func->SetRange(0.02, 8);
1782   func->DrawCopy("SAME");
1783   gPad->SetLogy();
1784
1785   // now fix exp. parameters and fit second part
1786   Double_t param0 = func->GetParameter(0);
1787   Double_t param1 = func->GetParameter(1);
1788   func->SetParameters(0, -1, 1, 1, 1);
1789   func->FixParameter(0, 0);
1790   func->FixParameter(1, -1);
1791   func->ReleaseParameter(2);
1792   func->ReleaseParameter(3);
1793   func->ReleaseParameter(4);
1794   func->SetParLimits(3, 1, 10);
1795   func->SetParLimits(4, 0, 10);
1796
1797   genePt->Fit(func, "0", "", 1.5, 4);
1798
1799   new TCanvas;
1800   genePt->DrawCopy("P");
1801   func->SetRange(0.02, 8);
1802   func->DrawCopy("SAME");
1803   gPad->SetLogy();
1804
1805   // fit both
1806   func->SetParameter(0, param0);
1807   func->SetParameter(1, param1);
1808   func->ReleaseParameter(0);
1809   func->ReleaseParameter(1);
1810
1811   new TCanvas;
1812   genePt->DrawCopy("P");
1813   func->SetRange(0.02, 8);
1814   func->DrawCopy("SAME");
1815   gPad->SetLogy();
1816
1817   genePt->Fit(func, "0", "", 0.2, 4);
1818
1819   TCanvas* canvas = new TCanvas("FitPt", "FitPt", 600, 400);
1820   canvas->SetGridx();
1821   canvas->SetGridy();
1822   canvas->SetRightMargin(0.05);
1823   canvas->SetTopMargin(0.05);
1824
1825   genePt->DrawCopy("P");
1826   func->SetRange(0.02, 8);
1827   func->DrawCopy("SAME");
1828   gPad->SetLogy();
1829
1830   TH1* first = (TH1*) func->GetHistogram()->Clone("first");
1831
1832   TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
1833
1834   TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
1835
1836   for (Int_t param=0; param<5; param++)
1837   {
1838     for (Int_t sign=0; sign<2; sign++)
1839     {
1840       TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign));  // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
1841       func2->SetParameters(func->GetParameters());
1842       //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
1843
1844       Float_t factor = ((sign == 0) ? 0.9 : 1.1);
1845       func2->SetParameter(param, func2->GetParameter(param) * factor);
1846       //func2->Print();
1847
1848       canvas->cd();
1849       func2->SetLineWidth(1);
1850       func2->SetLineColor(param + 2);
1851       func2->DrawCopy("SAME");
1852
1853       canvas2->cd();
1854       TH1* second = func2->GetHistogram();
1855       second->Divide(first);
1856       second->SetLineColor(param + 1);
1857       second->GetYaxis()->SetRangeUser(0, 2);
1858       second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
1859       second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
1860     }
1861   }
1862
1863   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1864
1865   file->Close();
1866 }
1867
1868 void SystematicpT()
1869 {
1870   gSystem->Load("libPWG0base");
1871
1872   TFile::Open("ptspectrum400.root");
1873   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1874   mult->LoadHistograms("Multiplicity");
1875
1876   AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
1877
1878   TH1* mcHist[12];
1879   TH1* result[12];
1880
1881   Int_t nParams = 1;
1882
1883   for (Int_t param=0; param<nParams; param++)
1884   {
1885     for (Int_t sign=0; sign<2; sign++)
1886     {
1887       // calculate result with systematic effect
1888       TFile* file = TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
1889       mult2->LoadHistograms("Multiplicity");
1890
1891       mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1892
1893       mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
1894       mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1895       //mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1896
1897       Int_t id = param * 2 + sign;
1898
1899       mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
1900       result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
1901
1902       TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
1903       DrawResultRatio(mcHist[id], result[id], tmp);
1904     }
1905   }
1906
1907   // calculate normal result
1908   TFile::Open("ptspectrum100_1.root");
1909   mult2->LoadHistograms("Multiplicity");
1910   TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc2");
1911
1912   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
1913
1914   mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
1915   //mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
1916
1917   TH1* result2 = mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
1918
1919   DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
1920
1921   DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
1922
1923   TCanvas* canvas = DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps");
1924
1925   TFile::Open("bayesianUncertainty_pt_400_100.root");
1926   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
1927   errorHist->SetLineColor(1);
1928   errorHist->SetLineWidth(2);
1929   TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
1930   for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
1931   {
1932     errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
1933     errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
1934   }
1935
1936   errorHist->DrawCopy("SAME");
1937   errorHist2->DrawCopy("SAME");
1938
1939   canvas->SaveAs(canvas->GetName());
1940
1941   DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
1942
1943   // does not make sense: mc is different
1944   //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
1945 }