PWG0/dNdEta/drawPlots.C

   1 /* $Id$ */
   2
   3 #if !defined(__CINT__) || defined(__MAKECINT__)
   4
   5 #include "AliPWG0Helper.h"
   6 #include "dNdEtaAnalysis.h"
   7 #include "AlidNdEtaCorrection.h"
   8
   9 #include <TCanvas.h>
  10 #include <TFile.h>
  11 #include <TH1.h>
  12 #include <TH2F.h>
  13 #include <TH3F.h>
  14 #include <TLine.h>
  15 #include <TSystem.h>
  16
  17 #endif
  18
  19 Int_t gMax = 5;
  20
  21 extern TSystem* gSystem;
  22
  23 void SetRanges(TAxis* axis)
  24 {
  25   if (strcmp(axis->GetTitle(), "#eta") == 0)
  26     axis->SetRangeUser(-1.7999, 1.7999);
  27   if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
  28     axis->SetRangeUser(0, 9.9999);
  29   if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
  30     axis->SetRangeUser(-15, 14.9999);
  31   if (strcmp(axis->GetTitle(), "Ntracks") == 0)
  32     axis->SetRangeUser(0, 99.9999);
  33 }
  34
  35 void SetRanges(TH1* hist)
  36 {
  37   SetRanges(hist->GetXaxis());
  38   SetRanges(hist->GetYaxis());
  39   SetRanges(hist->GetZaxis());
  40 }
  41
  42
  43 void Prepare3DPlot(TH3* hist)
  44 {
  45   hist->GetXaxis()->SetTitleOffset(1.5);
  46   hist->GetYaxis()->SetTitleOffset(1.5);
  47   hist->GetZaxis()->SetTitleOffset(1.5);
  48
  49   hist->SetStats(kFALSE);
  50 }
  51
  52 void Prepare2DPlot(TH2* hist)
  53 {
  54   hist->SetStats(kFALSE);
  55   hist->GetYaxis()->SetTitleOffset(1.4);
  56
  57   hist->SetMinimum(0);
  58   hist->SetMaximum(gMax);
  59
  60   SetRanges(hist);
  61 }
  62
  63 void Prepare1DPlot(TH1* hist)
  64 {
  65   hist->SetLineWidth(2);
  66   hist->SetStats(kFALSE);
  67
  68   hist->GetXaxis()->SetLabelOffset(0.02);
  69   hist->GetXaxis()->SetTitleOffset(1.3);
  70   hist->GetYaxis()->SetTitleOffset(1.3);
  71
  72   SetRanges(hist);
  73 }
  74
  75 void InitPad()
  76 {
  77   gPad->Range(0, 0, 1, 1);
  78   gPad->SetLeftMargin(0.15);
  79   //gPad->SetRightMargin(0.05);
  80   //gPad->SetTopMargin(0.13);
  81   gPad->SetBottomMargin(0.12);
  82
  83   gPad->SetGridx();
  84   gPad->SetGridy();
  85 }
  86
  87 void InitPadCOLZ()
  88 {
  89   gPad->Range(0, 0, 1, 1);
  90   gPad->SetRightMargin(0.15);
  91   gPad->SetLeftMargin(0.12);
  92
  93   gPad->SetGridx();
  94   gPad->SetGridy();
  95 }
  96
  97 void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
  98 {
  99   gSystem->Load("libPWG0base");
 100
 101   TFile::Open(esdFile);
 102   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 103   fdNdEtaAnalysisESD->LoadHistograms();
 104
 105   TFile::Open(mcFile);
 106   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
 107   fdNdEtaAnalysisMC->LoadHistograms();
 108   //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 109
 110   for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
 111     fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
 112
 113   fdNdEtaAnalysisESD->DrawHistograms();
 114 }
 115
 116 void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 117 {
 118   gSystem->Load("libPWG0base");
 119
 120   const char* ESDfolder = 0;
 121
 122   if (plot == 0) // all
 123     ESDfolder = "dndeta";
 124   else if (plot == 1) // mb
 125     ESDfolder = "dndeta_mb";
 126   else if (plot == 2) // mb vtx
 127     ESDfolder = "dndeta_mbvtx";
 128
 129   TFile::Open("analysis_esd.root");
 130   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 131   fdNdEtaAnalysisESD->LoadHistograms();
 132
 133   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 134   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 135
 136   TH2F* hist = 0;
 137
 138   if (plot == 0) // all
 139     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
 140   else if (plot == 1) // mb
 141     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
 142   else if (plot == 2) // mb vtx
 143     hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
 144
 145   TH1* proj = hist->ProjectionX();
 146
 147   TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
 148   for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
 149   {
 150     Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
 151     if (value != 0)
 152     {
 153       printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
 154       vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
 155     }
 156   }
 157
 158   new TCanvas;
 159   vertex->DrawCopy();
 160 }
 161
 162 void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 163 {
 164   gSystem->Load("libPWG0base");
 165
 166   TFile::Open("analysis_esd.root");
 167   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
 168   fdNdEtaAnalysisESD->LoadHistograms();
 169
 170   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 171   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 172
 173   //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
 174   //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 175
 176   TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
 177   TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 178
 179   TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
 180
 181   new TCanvas;
 182   histESD->Draw();
 183
 184   new TCanvas;
 185   histCorr->Draw();
 186
 187   for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
 188     for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
 189       for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
 190       {
 191         Float_t value1 = histESD->GetBinContent(x, y, z);
 192         Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
 193
 194         if (value2 > 0 && value1 > 0)
 195         {
 196           printf("%f %f %f\n", value1, value2, value1 / value2);
 197           diff->Fill(value1 / value2);
 198         }
 199       }
 200
 201   new TCanvas;
 202   diff->Draw();
 203 }
 204
 205 void dNdEta(Bool_t onlyESD = kFALSE)
 206 {
 207   TFile* file = TFile::Open("analysis_esd.root");
 208   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
 209   TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
 210   TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
 211   TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
 212   TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
 213   TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
 214   TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
 215
 216   Prepare1DPlot(histESD);
 217   Prepare1DPlot(histESDMB);
 218   Prepare1DPlot(histESDMBVtx);
 219
 220   Prepare1DPlot(histESDNoPt);
 221   Prepare1DPlot(histESDMBNoPt);
 222   Prepare1DPlot(histESDMBVtxNoPt);
 223   Prepare1DPlot(histESDMBTracksNoPt);
 224
 225   histESD->SetLineWidth(0);
 226   histESDMB->SetLineWidth(0);
 227   histESDMBVtx->SetLineWidth(0);
 228
 229   histESDNoPt->SetLineWidth(0);
 230   histESDMBNoPt->SetLineWidth(0);
 231   histESDMBVtxNoPt->SetLineWidth(0);
 232
 233   histESD->SetMarkerColor(1);
 234   histESDMB->SetMarkerColor(2);
 235   histESDMBVtx->SetMarkerColor(3);
 236
 237   histESDNoPt->SetMarkerColor(1);
 238   histESDMBNoPt->SetMarkerColor(2);
 239   histESDMBVtxNoPt->SetMarkerColor(3);
 240   histESDMBTracksNoPt->SetMarkerColor(4);
 241
 242   histESD->SetMarkerStyle(20);
 243   histESDMB->SetMarkerStyle(21);
 244   histESDMBVtx->SetMarkerStyle(22);
 245
 246   histESDNoPt->SetMarkerStyle(20);
 247   histESDMBNoPt->SetMarkerStyle(21);
 248   histESDMBVtxNoPt->SetMarkerStyle(22);
 249   histESDMBTracksNoPt->SetMarkerStyle(23);
 250
 251   TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 1000, 3.1, histESDMBVtx->GetMaximum() * 1.1);
 252   Prepare1DPlot(dummy);
 253   dummy->SetStats(kFALSE);
 254   dummy->SetXTitle("#eta");
 255   dummy->SetYTitle("dN_{ch}/d#eta");
 256   dummy->GetYaxis()->SetTitleOffset(1);
 257
 258   Float_t etaLimit = 0.7999;
 259
 260   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 261   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 262   histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 263
 264   histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 265   histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 266   histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 267   histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 268
 269   /*TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
 270
 271   dummy->DrawCopy();
 272   histESDMBVtx->Draw("SAME");
 273   histESDMB->Draw("SAME");
 274   histESD->Draw("SAME");
 275
 276   canvas->SaveAs("dNdEta1.gif");
 277   canvas->SaveAs("dNdEta1.eps");*/
 278
 279   if (onlyESD)
 280     return;
 281
 282   gSystem->Load("libPWG0base");
 283
 284   TFile* file2 = TFile::Open("analysis_mc.root");
 285   TH1* histMC = (TH1*) file2->Get("dndeta/dNdEta_corrected")->Clone("cloned");
 286   TH1* histMCTr = (TH1*) file2->Get("dndetaTr/dNdEta_corrected")->Clone("cloned2");
 287   TH1* histMCTrVtx = (TH1*) file2->Get("dndetaTrVtx/dNdEta_corrected")->Clone("cloned3");
 288
 289   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 290   fdNdEtaAnalysis->LoadHistograms();
 291   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 292   TH1* histMCPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 293
 294   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
 295   fdNdEtaAnalysis->LoadHistograms();
 296   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 297   TH1* histMCTrPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 298
 299   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
 300   fdNdEtaAnalysis->LoadHistograms();
 301   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 302   TH1* histMCTrVtxPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 303
 304   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
 305   fdNdEtaAnalysis->LoadHistograms();
 306   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 307   TH1* histMCTracksPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 308
 309   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
 310
 311   Prepare1DPlot(histMC);
 312   Prepare1DPlot(histMCTr);
 313   Prepare1DPlot(histMCTrVtx);
 314
 315   Prepare1DPlot(histMCPtCut);
 316   Prepare1DPlot(histMCTrPtCut);
 317   Prepare1DPlot(histMCTrVtxPtCut);
 318   if (histMCTracksPtCut)
 319     Prepare1DPlot(histMCTracksPtCut);
 320
 321   histMC->SetLineColor(1);
 322   histMCTr->SetLineColor(2);
 323   histMCTrVtx->SetLineColor(3);
 324
 325   histMCPtCut->SetLineColor(1);
 326   histMCTrPtCut->SetLineColor(2);
 327   histMCTrVtxPtCut->SetLineColor(3);
 328   if (histMCTracksPtCut)
 329     histMCTracksPtCut->SetLineColor(4);
 330
 331   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
 332   Prepare1DPlot(dummy2);
 333   dummy2->GetYaxis()->SetRangeUser(0, histESDMBVtx->GetMaximum() * 1.1);
 334
 335   dummy2->DrawCopy();
 336   histMC->Draw("SAME");
 337   histMCTr->Draw("SAME");
 338   histMCTrVtx->Draw("SAME");
 339   histESD->Draw("SAME");
 340   histESDMB->Draw("SAME");
 341   histESDMBVtx->Draw("SAME");
 342   histESDNoPt->Draw("SAME");
 343   histESDMBNoPt->Draw("SAME");
 344   histESDMBVtxNoPt->Draw("SAME");
 345   histESDMBTracksNoPt->Draw("SAME");
 346   histMCPtCut->Draw("SAME");
 347   histMCTrPtCut->Draw("SAME");
 348   histMCTrVtxPtCut->Draw("SAME");
 349   if (histMCTracksPtCut)
 350     histMCTracksPtCut->Draw("SAME");
 351
 352   canvas2->SaveAs("dNdEta2.gif");
 353   canvas2->SaveAs("dNdEta2.eps");
 354
 355   TH1* ratio = (TH1*) histMC->Clone("ratio");
 356   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
 357
 358   ratio->Divide(histESD);
 359   ratioNoPt->Divide(histESDNoPt);
 360
 361   ratio->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
 362
 363   ratio->SetLineColor(1);
 364   ratioNoPt->SetLineColor(2);
 365
 366   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 600);
 367   canvas3->Range(0, 0, 1, 1);
 368   //canvas3->Divide(1, 2, 0, 0);
 369
 370   //canvas3->cd(1);
 371   TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
 372   pad1->Draw();
 373
 374   TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
 375   pad2->Draw();
 376
 377   pad1->SetRightMargin(0.05);
 378   pad2->SetRightMargin(0.05);
 379
 380   // no border between them
 381   pad1->SetBottomMargin(0);
 382   pad2->SetTopMargin(0);
 383
 384   pad1->cd();
 385
 386   TLegend* legend = new TLegend(0.4, 0.05, 0.65, 0.3);
 387   legend->SetFillColor(0);
 388   legend->AddEntry(histESDMBVtx, "triggered, vertex");
 389   legend->AddEntry(histESDMB, "triggered");
 390   legend->AddEntry(histESD, "all events");
 391   legend->AddEntry(histMC, "MC prediction");
 392
 393   dummy->GetXaxis()->SetLabelSize(0.06);
 394   dummy->GetYaxis()->SetLabelSize(0.06);
 395   dummy->GetXaxis()->SetTitleSize(0.06);
 396   dummy->GetYaxis()->SetTitleSize(0.06);
 397   dummy->GetYaxis()->SetTitleOffset(0.7);
 398   dummy->DrawCopy();
 399   histESDMBVtx->Draw("SAME");
 400   histESDMB->Draw("SAME");
 401   histESD->Draw("SAME");
 402   histMC->Draw("SAME");
 403
 404   legend->Draw();
 405
 406   pad2->cd();
 407   pad2->SetBottomMargin(0.15);
 408
 409   Float_t min = TMath::Min(0.961, ratio->GetMinimum() * 0.95);
 410   Float_t max = TMath::Max(1.049, ratio->GetMaximum() * 1.05);
 411
 412   TH1F dummy3("dummy3", ";#eta;Ratio: MC / ESD", 1, -1.5, 1.5);
 413   dummy3.SetStats(kFALSE);
 414   dummy3.SetBinContent(1, 1);
 415   dummy3.GetYaxis()->SetRangeUser(min, max);
 416   dummy3.SetLineWidth(2);
 417   dummy3.GetXaxis()->SetLabelSize(0.06);
 418   dummy3.GetYaxis()->SetLabelSize(0.06);
 419   dummy3.GetXaxis()->SetTitleSize(0.06);
 420   dummy3.GetYaxis()->SetTitleSize(0.06);
 421   dummy3.GetYaxis()->SetTitleOffset(0.7);
 422   dummy3.DrawCopy();
 423
 424   ratio->Draw("SAME");
 425
 426   //pad2->Draw();
 427
 428   canvas3->Modified();
 429
 430   canvas3->SaveAs("dNdEta.gif");
 431   canvas3->SaveAs("dNdEta.eps");
 432
 433   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
 434
 435   ratio->Draw();
 436   ratioNoPt->Draw("SAME");
 437
 438   TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
 439   legend->SetFillColor(0);
 440   legend->AddEntry(ratio, "mc/esd");
 441   legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
 442   legend->Draw();
 443 }
 444
 445 void ptSpectrum()
 446 {
 447   TFile* file = TFile::Open("analysis_esd.root");
 448   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
 449
 450   TFile* file2 = TFile::Open("analysis_mc.root");
 451   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
 452
 453   TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
 454   InitPad();
 455   gPad->SetLogy();
 456
 457   Prepare1DPlot(histMC);
 458   Prepare1DPlot(histESD);
 459
 460   histESD->SetTitle("");
 461   histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 462   histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
 463
 464   histMC->SetLineColor(kBlue);
 465   histESD->SetLineColor(kRed);
 466
 467   histESD->GetYaxis()->SetTitleOffset(1.5);
 468   histESD->GetXaxis()->SetRangeUser(0, 4.9999);
 469
 470   histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
 471
 472   histESD->Draw();
 473   histMC->Draw("SAME");
 474
 475   canvas->SaveAs("ptSpectrum.gif");
 476   canvas->SaveAs("ptSpectrum.eps");
 477 }
 478
 479 void ptCutoff()
 480 {
 481   gSystem->Load("libPWG0base");
 482
 483   TFile::Open("correction_map.root");
 484   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 485   dNdEtaCorrection->LoadHistograms();
 486
 487   dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, -100, kTRUE);
 488
 489   TH1* hist = dynamic_cast<TH1*> (gROOT->FindObject("generated_pt")->Clone("ptcutoff"));
 490
 491   hist->GetXaxis()->SetRangeUser(0, 0.9999);
 492   hist->SetMinimum(0);
 493
 494   hist->SetTitle("Generated Particles");
 495   Prepare1DPlot(hist);
 496
 497   TCanvas* canvas = new TCanvas("ptCutoff", "ptCutoff", 700, 500);
 498   hist->DrawCopy();
 499
 500   TLine* line = new TLine(0.3, 0 - hist->GetMaximum() * 0, 0.3, hist->GetMaximum() * 1.1);
 501   line->SetLineWidth(3);
 502   line->SetLineColor(kRed);
 503   line->Draw();
 504
 505   canvas->SaveAs("ptCutoff.gif");
 506   canvas->SaveAs("ptCutoff.eps");
 507
 508   TH1F* factor = new TH1F("factor", ";#eta;correction factor", 20, -1, 1.000001);
 509   factor->SetLineWidth(2);
 510   for (Float_t eta = -0.95; eta<1; eta += 0.1)
 511     factor->Fill(eta, 1.0 / dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, eta, kFALSE));
 512
 513   TCanvas* canvas = new TCanvas("ptCutoff_factor", "ptCutoff_factor", 700, 500);
 514   InitPad();
 515
 516   Prepare1DPlot(factor);
 517   factor->GetYaxis()->SetRangeUser(1, 2);
 518   factor->GetYaxis()->SetTitleOffset(1);
 519   factor->Draw();
 520
 521   canvas->SaveAs("ptCutoff_factor.eps");
 522 }
 523
 524 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
 525 {
 526   gSystem->Load("libPWG0base");
 527
 528   TFile::Open(fileName);
 529   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 530   dNdEtaCorrection->LoadHistograms();
 531
 532   TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
 533   TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
 534
 535   Prepare2DPlot(corrTrigger);
 536   corrTrigger->SetTitle("b) Trigger bias correction");
 537
 538   Prepare2DPlot(corrVtx);
 539   corrVtx->SetTitle("a) Vertex reconstruction correction");
 540
 541   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
 542   corrVtx->GetYaxis()->SetTitle("Multiplicity");
 543
 544   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
 545   canvas->Divide(2, 1);
 546
 547   canvas->cd(1);
 548   InitPadCOLZ();
 549   corrVtx->DrawCopy("COLZ");
 550
 551   canvas->cd(2);
 552   InitPadCOLZ();
 553   corrTrigger->DrawCopy("COLZ");
 554
 555   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
 556   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
 557
 558   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
 559   canvas->Divide(2, 1);
 560
 561   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
 562   corrVtx->GetYaxis()->SetRangeUser(0, 5);
 563
 564   canvas->cd(1);
 565   InitPadCOLZ();
 566   corrVtx->DrawCopy("COLZ");
 567
 568   canvas->cd(2);
 569   InitPadCOLZ();
 570   corrTrigger->DrawCopy("COLZ");
 571
 572   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
 573   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
 574 }
 575
 576 void TriggerBias(const char* fileName = "correction_map.root")
 577 {
 578   TFile* file = TFile::Open(fileName);
 579
 580   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
 581
 582   Prepare2DPlot(corr);
 583   corr->SetTitle("Trigger bias correction");
 584
 585   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
 586   InitPadCOLZ();
 587   corr->DrawCopy("COLZ");
 588
 589   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
 590   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
 591
 592   corr->GetYaxis()->SetRangeUser(0, 5);
 593
 594   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
 595   InitPadCOLZ();
 596   corr->DrawCopy("COLZ");
 597
 598   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
 599   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
 600 }
 601
 602 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
 603 {
 604   gSystem->Load("libPWG0base");
 605
 606   TFile* file = TFile::Open(fileName);
 607   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
 608   dNdEtaCorrection->LoadHistograms();
 609
 610   TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
 611   TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
 612
 613   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
 614   canvas->Divide(2, 1);
 615
 616   canvas->cd(1);
 617   InitPad();
 618
 619   Prepare1DPlot(hist);
 620   hist->SetTitle("");
 621   hist->GetYaxis()->SetTitle("correction factor");
 622   hist->GetYaxis()->SetRangeUser(1, 1.5);
 623   hist->GetYaxis()->SetTitleOffset(1.6);
 624   hist->Draw();
 625
 626   canvas->cd(2);
 627   InitPad();
 628
 629   Prepare1DPlot(hist2);
 630   hist2->SetTitle("");
 631   hist2->GetYaxis()->SetTitle("correction factor");
 632   hist2->GetXaxis()->SetRangeUser(0, 5);
 633   hist2->GetYaxis()->SetTitleOffset(1.6);
 634   hist2->GetXaxis()->SetTitle("multiplicity");
 635   hist2->Draw();
 636
 637   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
 638   pave->SetFillColor(0);
 639   pave->AddText("|z| < 10 cm");
 640   pave->Draw();
 641
 642   canvas->SaveAs("TriggerBias1D.eps");
 643 }
 644
 645 void VtxRecon()
 646 {
 647   TFile* file = TFile::Open("correction_map.root");
 648
 649   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
 650
 651   Prepare2DPlot(corr);
 652   corr->SetTitle("Vertex reconstruction correction");
 653
 654   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
 655   InitPadCOLZ();
 656   corr->DrawCopy("COLZ");
 657
 658   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
 659   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
 660
 661   corr->GetYaxis()->SetRangeUser(0, 5);
 662
 663   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
 664   InitPadCOLZ();
 665   corr->DrawCopy("COLZ");
 666
 667   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
 668   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
 669 }
 670
 671 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
 672 {
 673   gSystem->Load("libPWG0base");
 674
 675   TFile* file = TFile::Open(fileName);
 676   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
 677   dNdEtaCorrection->LoadHistograms();
 678
 679   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
 680   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
 681
 682   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
 683   canvas->Divide(2, 1);
 684
 685   canvas->cd(1);
 686   InitPad();
 687
 688   Prepare1DPlot(hist);
 689   hist->SetTitle("");
 690   hist->GetYaxis()->SetTitle("correction factor");
 691   hist->GetYaxis()->SetRangeUser(1, 1.8);
 692   hist->GetYaxis()->SetTitleOffset(1.6);
 693   hist->DrawCopy();
 694
 695   canvas->cd(2);
 696   InitPad();
 697
 698   Prepare1DPlot(hist2);
 699   hist2->SetTitle("");
 700   hist2->GetYaxis()->SetTitle("correction factor");
 701   hist2->GetXaxis()->SetRangeUser(0, 20);
 702   hist2->GetYaxis()->SetTitleOffset(1.6);
 703   hist2->GetXaxis()->SetTitle("multiplicity");
 704   hist2->Draw();
 705
 706   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
 707   pave->SetFillColor(0);
 708   pave->AddText("|z| < 10 cm");
 709   pave->Draw();
 710
 711   canvas->SaveAs("VtxRecon1D.eps");
 712
 713   Correction1DCreatePlots(fileName, folderName, 9.9, 2);
 714
 715   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
 716   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
 717
 718   Prepare1DPlot(corrX);
 719   Prepare1DPlot(corrZ);
 720
 721   corrX->GetYaxis()->SetTitleOffset(1.5);
 722   corrZ->GetYaxis()->SetTitleOffset(1.5);
 723
 724   corrX->SetTitle("a) z projection");
 725   corrZ->SetTitle("b) p_{T} projection");
 726
 727   corrX->GetYaxis()->SetTitle("correction factor");
 728   corrZ->GetYaxis()->SetTitle("correction factor");
 729
 730   corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
 731
 732   TString canvasName;
 733   canvasName.Form("VtxRecon1D_Track");
 734   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
 735   canvas->Divide(2, 1);
 736
 737   canvas->cd(1);
 738   InitPad();
 739   corrX->DrawCopy();
 740
 741   canvas->cd(2);
 742   InitPad();
 743   gPad->SetLogx();
 744   corrZ->Draw();
 745
 746   canvas->SaveAs("VtxRecon1D_Track.eps");
 747   canvas->SaveAs("VtxRecon1D_Track.gif");
 748 }
 749
 750 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
 751 {
 752   gSystem->Load("libPWG0base");
 753
 754   TFile::Open(fileName);
 755   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 756   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
 757
 758   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 759   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 760
 761   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 762   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 763
 764   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 765   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
 766   gene->GetXaxis()->SetRangeUser(-10, 10);
 767   meas->GetXaxis()->SetRangeUser(-10, 10);
 768
 769   Float_t eff1 = gene->Integral() / meas->Integral();
 770   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
 771
 772   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
 773
 774   gene->GetZaxis()->SetRangeUser(0.3, 10);
 775   meas->GetZaxis()->SetRangeUser(0.3, 10);
 776
 777   Float_t eff2 = gene->Integral() / meas->Integral();
 778   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
 779
 780   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
 781
 782   gene->GetZaxis()->SetRangeUser(0.3, 1);
 783   meas->GetZaxis()->SetRangeUser(0.3, 1);
 784
 785   Float_t eff3 = gene->Integral() / meas->Integral();
 786   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
 787
 788   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
 789 }
 790
 791 void Correction1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType = 0)
 792 {
 793   TFile::Open(fileName);
 794   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
 795   dNdEtaCorrection->LoadHistograms();
 796
 797   TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
 798   TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetMeasuredHistogram();
 799
 800   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
 801   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
 802   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 803   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
 804   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kTRUE);
 805   gene->GetYaxis()->SetRange(0, 0);
 806   meas->GetYaxis()->SetRange(0, 0);
 807
 808   gene->GetXaxis()->SetRangeUser(-10, 10);
 809   meas->GetXaxis()->SetRangeUser(-10, 10);
 810   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kTRUE);
 811   gene->GetZaxis()->SetRange(0, 0);
 812   meas->GetZaxis()->SetRange(0, 0);
 813
 814   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 815   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
 816   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kTRUE);
 817 }
 818
 819 void Correction1D(Int_t correctionType = 0, const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
 820 {
 821   gSystem->Load("libPWG0base");
 822
 823   Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType);
 824
 825   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
 826   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
 827   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
 828
 829   Prepare1DPlot(corrX);
 830   Prepare1DPlot(corrY);
 831   Prepare1DPlot(corrZ);
 832
 833   corrX->SetTitle("a) z projection");
 834   corrY->SetTitle("b) #eta projection");
 835   corrZ->SetTitle("c) p_{T} projection");
 836
 837   corrX->GetYaxis()->SetTitle("correction factor");
 838   corrY->GetYaxis()->SetTitle("correction factor");
 839   corrZ->GetYaxis()->SetTitle("correction factor");
 840
 841   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
 842
 843   TString canvasName;
 844   canvasName.Form("Correction1D_%s", folder);
 845   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
 846   canvas->Divide(3, 1);
 847
 848   canvas->cd(1);
 849   InitPad();
 850   corrX->DrawCopy();
 851
 852   canvas->cd(2);
 853   InitPad();
 854   corrY->Draw();
 855
 856   canvas->cd(3);
 857   InitPad();
 858   corrZ->Draw();
 859
 860   canvas->SaveAs(Form("Correction1D_%d_%s_%f.gif", correctionType, fileName, upperPtLimit));
 861   canvas->SaveAs(Form("Correction1D_%d_%s_%f.eps", correctionType, fileName, upperPtLimit));
 862 }
 863
 864 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
 865 {
 866   gSystem->Load("libPWG0base");
 867
 868   Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
 869
 870   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
 871   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
 872   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
 873
 874   Prepare1DPlot(corrX);
 875   Prepare1DPlot(corrY);
 876   Prepare1DPlot(corrZ);
 877
 878   corrX->SetTitle("a) z projection");
 879   corrY->SetTitle("a) #eta projection");
 880   corrZ->SetTitle("b) p_{T} projection");
 881
 882   corrY->GetYaxis()->SetTitle("correction factor");
 883   corrZ->GetYaxis()->SetTitle("correction factor");
 884
 885   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
 886
 887   TString canvasName;
 888   canvasName.Form("Track2Particle1D_%s", folder);
 889   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
 890   canvas->Divide(3, 1);
 891
 892   canvas->cd(1);
 893   InitPad();
 894   corrX->DrawCopy();
 895
 896   canvas->cd(2);
 897   InitPad();
 898   corrY->Draw();
 899
 900   canvas->cd(3);
 901   InitPad();
 902   corrZ->Draw();
 903
 904   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
 905   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
 906
 907   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
 908
 909   canvasName.Form("Track2Particle1D_%s_etapt", folder);
 910   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
 911   canvas->Divide(2, 1);
 912
 913   canvas->cd(1);
 914   InitPad();
 915   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
 916   corrY->GetYaxis()->SetRangeUser(1, 1.5);
 917   corrY->GetYaxis()->SetTitleOffset(1.5);
 918   corrY->DrawCopy();
 919   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
 920   pave->AddText("|z| < 10 cm");
 921   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
 922   pave->Draw();
 923
 924   canvas->cd(2);
 925   InitPad();
 926   gPad->SetLogx();
 927   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
 928   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
 929   corrZ->GetYaxis()->SetTitleOffset(1.5);
 930   corrZ->DrawCopy();
 931   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
 932   pave->AddText("|z| < 10 cm");
 933   pave->AddText("|#eta| < 0.8");
 934   pave->Draw();
 935
 936   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
 937   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
 938 }
 939
 940 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
 941 {
 942   gSystem->Load("libPWG0base");
 943
 944   // particle type
 945   for (Int_t particle=0; particle<4; ++particle)
 946   {
 947     TString dirName;
 948     dirName.Form("correction_%d", particle);
 949     Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
 950
 951     TString tmpx, tmpy, tmpz;
 952     tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
 953     tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
 954     tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
 955
 956     TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
 957     TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
 958     TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
 959
 960     Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
 961
 962     TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
 963     TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
 964     TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
 965
 966     posX->Divide(negX);
 967     posY->Divide(negY);
 968     posZ->Divide(negZ);
 969
 970     Prepare1DPlot(posX);
 971     Prepare1DPlot(posY);
 972     Prepare1DPlot(posZ);
 973
 974     Float_t min = 0.8;
 975     Float_t max = 1.2;
 976
 977     posX->SetMinimum(min);
 978     posX->SetMaximum(max);
 979     posY->SetMinimum(min);
 980     posY->SetMaximum(max);
 981     posZ->SetMinimum(min);
 982     posZ->SetMaximum(max);
 983
 984     posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
 985
 986     posX->GetYaxis()->SetTitleOffset(1.7);
 987     posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
 988     posY->GetYaxis()->SetTitleOffset(1.7);
 989     posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
 990     posZ->GetYaxis()->SetTitleOffset(1.7);
 991     posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
 992
 993     posZ->GetXaxis()->SetRangeUser(0, 1);
 994
 995     TString canvasName;
 996     canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
 997
 998     TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
 999     canvas->Divide(3, 1);
1000
1001     canvas->cd(1);
1002     InitPad();
1003     posX->DrawCopy();
1004
1005     canvas->cd(2);
1006     InitPad();
1007     posY->DrawCopy();
1008
1009     canvas->cd(3);
1010     InitPad();
1011     posZ->DrawCopy();
1012
1013     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1014     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1015   }
1016 }
1017
1018 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
1019 {
1020   TFile::Open(fileName);
1021   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1022   dNdEtaCorrection->LoadHistograms();
1023
1024   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1025   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1026
1027   gene->GetZaxis()->SetRangeUser(0.3, 10);
1028   meas->GetZaxis()->SetRangeUser(0.3, 10);
1029   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
1030   gene->GetZaxis()->SetRange(0, 0);
1031   meas->GetZaxis()->SetRange(0, 0);
1032
1033   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1034   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1035   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
1036   gene->GetYaxis()->SetRange(0, 0);
1037   meas->GetYaxis()->SetRange(0, 0);
1038
1039   gene->GetXaxis()->SetRangeUser(-10, 10);
1040   meas->GetXaxis()->SetRangeUser(-10, 10);
1041   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
1042   gene->GetXaxis()->SetRange(0, 0);
1043   meas->GetXaxis()->SetRange(0, 0);
1044 }
1045
1046 void Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1047 {
1048   gSystem->Load("libPWG0base");
1049
1050   Track2Particle2DCreatePlots(fileName);
1051
1052   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1053   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1054   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
1055
1056   /* this reads them from the file
1057   TH2* corrYX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx"));
1058   TH2* corrZX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx"));
1059   TH2* corrZY = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy"));*/
1060
1061   Prepare2DPlot(corrYX);
1062   Prepare2DPlot(corrZX);
1063   Prepare2DPlot(corrZY);
1064
1065   const char* title = "";
1066   corrYX->SetTitle(title);
1067   corrZX->SetTitle(title);
1068   corrZY->SetTitle(title);
1069
1070   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1071   canvas->Divide(3, 1);
1072
1073   canvas->cd(1);
1074   InitPadCOLZ();
1075   corrYX->Draw("COLZ");
1076
1077   canvas->cd(2);
1078   InitPadCOLZ();
1079   corrZX->Draw("COLZ");
1080
1081   canvas->cd(3);
1082   InitPadCOLZ();
1083   corrZY->Draw("COLZ");
1084
1085   canvas->SaveAs(Form("Track2Particle2D_%d.gif", gMax));
1086   canvas->SaveAs(Form("Track2Particle2D_%d.eps", gMax));
1087 }
1088
1089 void CompareTrack2Particle2D()
1090 {
1091   gSystem->Load("libPWG0base");
1092
1093   Track2Particle2DCreatePlots("correction_maponly-positive.root");
1094
1095   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("pos_yx"));
1096   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("pos_zx"));
1097   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("pos_zy"));
1098
1099   Track2Particle2DCreatePlots("correction_maponly-negative.root");
1100
1101   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("neg_yx"));
1102   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("neg_zx"));
1103   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("neg_zy"));
1104
1105   posYX->Divide(negYX);
1106   posZX->Divide(negZX);
1107   posZY->Divide(negZY);
1108
1109   Prepare2DPlot(posYX);
1110   Prepare2DPlot(posZX);
1111   Prepare2DPlot(posZY);
1112
1113   Float_t min = 0.8;
1114   Float_t max = 1.2;
1115
1116   posYX->SetMinimum(min);
1117   posYX->SetMaximum(max);
1118   posZX->SetMinimum(min);
1119   posZX->SetMaximum(max);
1120   posZY->SetMinimum(min);
1121   posZY->SetMaximum(max);
1122
1123   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1124   canvas->Divide(3, 1);
1125
1126   canvas->cd(1);
1127   InitPadCOLZ();
1128   posYX->Draw("COLZ");
1129
1130   canvas->cd(2);
1131   InitPadCOLZ();
1132   posZX->Draw("COLZ");
1133
1134   canvas->cd(3);
1135   InitPadCOLZ();
1136   posZY->Draw("COLZ");
1137
1138   canvas->SaveAs("CompareTrack2Particle2D.gif");
1139   canvas->SaveAs("CompareTrack2Particle2D.eps");
1140 }
1141
1142 void Track2Particle3D()
1143 {
1144   // get left margin proper
1145
1146   TFile* file = TFile::Open("correction_map.root");
1147
1148   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1149
1150   corr->SetTitle("Correction Factor");
1151   SetRanges(corr->GetZaxis());
1152
1153   Prepare3DPlot(corr);
1154
1155   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1156   canvas->SetTheta(29.428);
1157   canvas->SetPhi(16.5726);
1158
1159   corr->Draw();
1160
1161   canvas->SaveAs("Track2Particle3D.gif");
1162   canvas->SaveAs("Track2Particle3D.eps");
1163 }
1164
1165 void Track2Particle3DAll()
1166 {
1167   TFile* file = TFile::Open("correction_map.root");
1168
1169   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1170   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1171   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1172
1173   gene->SetTitle("Generated Particles");
1174   meas->SetTitle("Measured Tracks");
1175   corr->SetTitle("Correction Factor");
1176
1177   Prepare3DPlot(gene);
1178   Prepare3DPlot(meas);
1179   Prepare3DPlot(corr);
1180
1181   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1182   canvas->Divide(3, 1);
1183
1184   canvas->cd(1);
1185   InitPad();
1186   gene->Draw();
1187
1188   canvas->cd(2);
1189   meas->Draw();
1190
1191   canvas->cd(3);
1192   corr->Draw();
1193
1194   canvas->SaveAs("Track2Particle3DAll.gif");
1195   canvas->SaveAs("Track2Particle3DAll.eps");
1196 }
1197
1198 void MultiplicityMC(Int_t xRangeMax = 50)
1199 {
1200   TFile* file = TFile::Open("multiplicityMC.root");
1201
1202   if (!file)
1203   {
1204     printf("multiplicityMC.root could not be opened.\n");
1205     return;
1206   }
1207
1208   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1209   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1210   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1211
1212   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1213   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1214   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1215   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1216   {
1217     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1218     proj->Fit("gaus", "0");
1219     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1220     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1221
1222     continue;
1223
1224     // draw for debugging
1225     new TCanvas;
1226     proj->DrawCopy();
1227     proj->GetFunction("gaus")->DrawCopy("SAME");
1228   }
1229
1230   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1231
1232   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1233   {
1234     Float_t mean = correction->GetBinContent(i);
1235     Float_t width = correctionWidth->GetBinContent(i);
1236
1237     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1238     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1239     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1240
1241     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1242     {
1243       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1244     }
1245   }
1246
1247   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1248   fMultiplicityESDCorrectedRebinned->Rebin(10);
1249   fMultiplicityESDCorrectedRebinned->Scale(0.1);
1250
1251   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1252   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1253   ratio->Divide(fMultiplicityMC);
1254
1255   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1256   ratio2->Divide(fMultiplicityMC);
1257
1258   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1259   canvas->Divide(3, 2);
1260
1261   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1262   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1263   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1264   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1265   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1266   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1267   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1268
1269   canvas->cd(1); //InitPad();
1270   fMultiplicityESD->Draw();
1271   fMultiplicityMC->SetLineColor(2);
1272   fMultiplicityMC->Draw("SAME");
1273
1274   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1275   legend->AddEntry(fMultiplicityESD, "ESD");
1276   legend->AddEntry(fMultiplicityMC, "MC");
1277   legend->Draw();
1278
1279   canvas->cd(2);
1280   fCorrelation->Draw("COLZ");
1281
1282   canvas->cd(3);
1283   correction->Draw();
1284   //correction->Fit("pol1");
1285   correctionWidth->SetLineColor(2);
1286   correctionWidth->Draw("SAME");
1287
1288   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1289   legend->AddEntry(correction, "#bar{x}");
1290   legend->AddEntry(correctionWidth, "#sigma");
1291   legend->Draw();
1292
1293   canvas->cd(4);
1294   ratio->Draw();
1295
1296   ratio2->SetLineColor(2);
1297   ratio2->Draw("SAME");
1298
1299   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1300   legend->AddEntry(ratio, "uncorrected");
1301   legend->AddEntry(ratio2, "corrected");
1302   legend->Draw();
1303
1304   canvas->cd(5);
1305   fMultiplicityESDCorrected->SetLineColor(kBlue);
1306   fMultiplicityESDCorrected->Draw();
1307   fMultiplicityMC->Draw("SAME");
1308   fMultiplicityESD->Draw("SAME");
1309
1310   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1311   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1312   legend->AddEntry(fMultiplicityMC, "MC");
1313   legend->AddEntry(fMultiplicityESD, "ESD");
1314   legend->Draw();
1315
1316   canvas->cd(6);
1317   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
1318   fMultiplicityESDCorrectedRebinned->Draw();
1319   fMultiplicityMC->Draw("SAME");
1320
1321   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1322   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1323   legend->AddEntry(fMultiplicityMC, "MC");
1324   legend->Draw();
1325
1326   canvas->SaveAs("MultiplicityMC.gif");
1327 }
1328
1329 void MultiplicityESD()
1330 {
1331   TFile* file = TFile::Open("multiplicityESD.root");
1332
1333   if (!file)
1334   {
1335     printf("multiplicityESD.root could not be opened.\n");
1336     return;
1337   }
1338
1339   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1340
1341   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1342
1343   fMultiplicityESD->Draw();
1344 }
1345
1346 void drawPlots(Int_t max)
1347 {
1348   gMax = max;
1349
1350   ptCutoff();
1351   TriggerBias();
1352   VtxRecon();
1353   Track2Particle2D();
1354   Track2Particle3D();
1355   ptSpectrum();
1356   dNdEta();
1357 }
1358
1359 void drawPlots()
1360 {
1361   drawPlots(5);
1362   drawPlots(2);
1363 }