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);
 109
 110   for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
 111     fdNdEtaAnalysisESD->GetdNdEtaHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaHistogram(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
 215   TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
 216
 217   Prepare1DPlot(histESD);
 218   Prepare1DPlot(histESDMB);
 219   Prepare1DPlot(histESDMBVtx);
 220
 221   Prepare1DPlot(histESDNoPt);
 222   Prepare1DPlot(histESDMBNoPt);
 223   Prepare1DPlot(histESDMBVtxNoPt);
 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(kRed);
 234   histESDMB->SetMarkerColor(kBlue);
 235   histESDMBVtx->SetMarkerColor(103);
 236
 237   histESDNoPt->SetMarkerColor(kRed);
 238   histESDMBNoPt->SetMarkerColor(kBlue);
 239   histESDMBVtxNoPt->SetMarkerColor(103);
 240
 241   histESD->SetMarkerStyle(20);
 242   histESDMB->SetMarkerStyle(21);
 243   histESDMBVtx->SetMarkerStyle(22);
 244
 245   histESDNoPt->SetMarkerStyle(20);
 246   histESDMBNoPt->SetMarkerStyle(21);
 247   histESDMBVtxNoPt->SetMarkerStyle(22);
 248
 249   TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 1000, 0.1, histESDMBVtx->GetMaximum() * 1.1);
 250   Prepare1DPlot(dummy);
 251   dummy->SetStats(kFALSE);
 252   dummy->SetXTitle("#eta");
 253   dummy->SetYTitle("dN_{ch}/d#eta");
 254   dummy->GetYaxis()->SetTitleOffset(1);
 255
 256   Float_t etaLimit = 0.7999;
 257
 258   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 259   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 260   histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 261
 262   histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 263   histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 264   histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 265
 266   dummy->DrawCopy();
 267   histESDMBVtx->Draw("SAME");
 268   histESDMB->Draw("SAME");
 269   histESD->Draw("SAME");
 270
 271   canvas->SaveAs("dNdEta1.gif");
 272   canvas->SaveAs("dNdEta1.eps");
 273
 274   if (onlyESD)
 275     return;
 276
 277   gSystem->Load("libPWG0base");
 278
 279   TFile* file2 = TFile::Open("analysis_mc.root");
 280   TH1* histMC = (TH1*) file2->Get("dndeta/dNdEta_corrected")->Clone("cloned");
 281   TH1* histMCTr = (TH1*) file2->Get("dndetaTr/dNdEta_corrected")->Clone("cloned2");
 282   TH1* histMCTrVtx = (TH1*) file2->Get("dndetaTrVtx/dNdEta_corrected")->Clone("cloned3");
 283
 284   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 285   fdNdEtaAnalysis->LoadHistograms();
 286   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 287   TH1* histMCPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 288
 289   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
 290   fdNdEtaAnalysis->LoadHistograms();
 291   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 292   TH1* histMCTrPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 293
 294   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
 295   fdNdEtaAnalysis->LoadHistograms();
 296   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 297   TH1* histMCTrVtxPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 298
 299   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
 300
 301   Prepare1DPlot(histMC);
 302   Prepare1DPlot(histMCTr);
 303   Prepare1DPlot(histMCTrVtx);
 304
 305   Prepare1DPlot(histMCPtCut);
 306   Prepare1DPlot(histMCTrPtCut);
 307   Prepare1DPlot(histMCTrVtxPtCut);
 308
 309   histMC->SetLineColor(kBlue);
 310   histMCTr->SetLineColor(kBlue);
 311   histMCTrVtx->SetLineColor(kBlue);
 312
 313   histMCPtCut->SetLineColor(kBlue);
 314   histMCTrPtCut->SetLineColor(kBlue);
 315   histMCTrVtxPtCut->SetLineColor(kBlue);
 316
 317   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
 318   Prepare1DPlot(dummy2);
 319   dummy2->GetYaxis()->SetRangeUser(0, histESDMBVtx->GetMaximum() * 1.1);
 320
 321   dummy2->DrawCopy();
 322   histMC->Draw("SAME");
 323   histMCTr->Draw("SAME");
 324   histMCTrVtx->Draw("SAME");
 325   histESD->Draw("SAME");
 326   histESDMB->Draw("SAME");
 327   histESDMBVtx->Draw("SAME");
 328   histESDNoPt->Draw("SAME");
 329   histESDMBNoPt->Draw("SAME");
 330   histESDMBVtxNoPt->Draw("SAME");
 331   histMCPtCut->Draw("SAME");
 332   histMCTrPtCut->Draw("SAME");
 333   histMCTrVtxPtCut->Draw("SAME");
 334
 335   canvas2->SaveAs("dNdEta2.gif");
 336   canvas2->SaveAs("dNdEta2.eps");
 337
 338   TH1* ratio = (TH1*) histMC->Clone("ratio");
 339   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
 340
 341   ratio->Divide(histESD);
 342   ratioNoPt->Divide(histESDNoPt);
 343
 344   ratio->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
 345
 346   ratio->SetLineColor(1);
 347   ratioNoPt->SetLineColor(2);
 348
 349   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 700);
 350   canvas3->Range(0, 0, 1, 1);
 351   //canvas3->Divide(1, 2, 0, 0);
 352
 353   //canvas3->cd(1);
 354   TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.4, 0.98, 0.98);
 355   pad1->Draw();
 356
 357   TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.4);
 358   pad2->Draw();
 359
 360   pad1->SetRightMargin(0.05);
 361   pad2->SetRightMargin(0.05);
 362
 363   // no border between them
 364   pad1->SetBottomMargin(0);
 365   pad2->SetTopMargin(0);
 366
 367   pad1->cd();
 368
 369   TLegend* legend = new TLegend(0.4, 0.2, 0.65, 0.4);
 370   legend->SetFillColor(0);
 371   legend->AddEntry(histESDMBVtx, "triggered, vertex");
 372   legend->AddEntry(histESDMB, "triggered");
 373   legend->AddEntry(histESD, "all events");
 374   legend->AddEntry(histMC, "MC prediction");
 375
 376   dummy->DrawCopy();
 377   histESDMBVtx->Draw("SAME");
 378   histESDMB->Draw("SAME");
 379   histESD->Draw("SAME");
 380   histMC->Draw("SAME");
 381
 382   legend->Draw();
 383
 384   pad2->cd();
 385   pad2->SetBottomMargin(0.15);
 386
 387   TH1F dummy3("dummy3", ";#eta;Ratio: MC / ESD", 1, -1.5, 1.5);
 388   dummy3.SetStats(kFALSE);
 389   dummy3.SetBinContent(1, 1);
 390   dummy3.GetYaxis()->SetRangeUser(0.961, 1.049);
 391   dummy3.SetLineWidth(2);
 392   dummy3.GetXaxis()->SetLabelSize(0.06);
 393   dummy3.GetYaxis()->SetLabelSize(0.06);
 394   dummy3.GetXaxis()->SetTitleSize(0.06);
 395   dummy3.GetYaxis()->SetTitleSize(0.06);
 396   dummy3.GetYaxis()->SetTitleOffset(0.7);
 397   dummy3.DrawCopy();
 398
 399   ratio->Draw("SAME");
 400
 401   //pad2->Draw();
 402
 403   canvas3->Modified();
 404
 405   canvas3->SaveAs("dNdEta.gif");
 406   canvas3->SaveAs("dNdEta.eps");
 407
 408   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
 409
 410   ratio->Draw();
 411   ratioNoPt->Draw("SAME");
 412
 413   TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
 414   legend->SetFillColor(0);
 415   legend->AddEntry(ratio, "mc/esd");
 416   legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
 417   legend->Draw();
 418 }
 419
 420 void ptSpectrum()
 421 {
 422   TFile* file = TFile::Open("analysis_esd.root");
 423   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
 424
 425   TFile* file2 = TFile::Open("analysis_mc.root");
 426   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
 427
 428   TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
 429   InitPad();
 430   gPad->SetLogy();
 431
 432   Prepare1DPlot(histMC);
 433   Prepare1DPlot(histESD);
 434
 435   histESD->SetTitle("");
 436   histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 437   histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
 438
 439   histMC->SetLineColor(kBlue);
 440   histESD->SetLineColor(kRed);
 441
 442   histESD->GetYaxis()->SetTitleOffset(1.5);
 443   histESD->GetXaxis()->SetRangeUser(0, 4.9999);
 444
 445   histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
 446
 447   histESD->Draw();
 448   histMC->Draw("SAME");
 449
 450   canvas->SaveAs("ptSpectrum.gif");
 451   canvas->SaveAs("ptSpectrum.eps");
 452 }
 453
 454 void ptCutoff()
 455 {
 456   gSystem->Load("libPWG0base");
 457
 458   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 459   dNdEtaCorrection->LoadHistograms("correction_map.root","dndeta_correction");
 460
 461   dNdEtaCorrection->GetMeasuredFraction(0.3, -100, kTRUE);
 462
 463   TH1* hist = dynamic_cast<TH1*> (gROOT->FindObject("gene_dndeta_correction_nTrackToNPart_pt")->Clone("ptcutoff"));
 464
 465   hist->GetXaxis()->SetRangeUser(0, 0.9999);
 466   hist->SetMinimum(0);
 467
 468   hist->SetTitle("Generated Particles");
 469   Prepare1DPlot(hist);
 470
 471   TCanvas* canvas = new TCanvas("ptCutoff", "ptCutoff", 700, 500);
 472   hist->DrawCopy();
 473
 474   TLine* line = new TLine(0.3, 0 - hist->GetMaximum() * 0, 0.3, hist->GetMaximum() * 1.1);
 475   line->SetLineWidth(3);
 476   line->SetLineColor(kRed);
 477   line->Draw();
 478
 479   canvas->SaveAs("ptCutoff.gif");
 480   canvas->SaveAs("ptCutoff.eps");
 481
 482   TH1F* factor = new TH1F("factor", ";#eta;correction factor", 10, -1, 1.000001);
 483   for (Float_t eta = -0.9; eta<1; eta += 0.2)
 484     factor->Fill(eta, 1.0 / dNdEtaCorrection->GetMeasuredFraction(0.3, eta, kFALSE));
 485
 486   TCanvas* canvas = new TCanvas("ptCutoff_factor", "ptCutoff_factor", 700, 500);
 487   InitPad();
 488
 489   Prepare1DPlot(factor);
 490   factor->GetYaxis()->SetRangeUser(1, 2);
 491   factor->GetYaxis()->SetTitleOffset(1);
 492   factor->Draw();
 493
 494   canvas->SaveAs("ptCutoff_factor.eps");
 495 }
 496
 497 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
 498 {
 499   TFile* file = TFile::Open(fileName);
 500
 501   TH2* corrTrigger = dynamic_cast<TH2*> (file->Get(Form("%s/corr_%s_trigger", folder, folder)));
 502   TH2* corrVtx = dynamic_cast<TH2*> (file->Get(Form("%s/corr_%s_vtxReco", folder, folder)));
 503
 504   Prepare2DPlot(corrTrigger);
 505   corrTrigger->SetTitle("a) Trigger bias correction");
 506
 507   Prepare2DPlot(corrVtx);
 508   corrVtx->SetTitle("b) Vertex reconstruction correction");
 509
 510   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
 511   corrVtx->GetYaxis()->SetTitle("Multiplicity");
 512
 513   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
 514   canvas->Divide(2, 1);
 515
 516   canvas->cd(1);
 517   InitPadCOLZ();
 518   corrTrigger->DrawCopy("COLZ");
 519
 520   canvas->cd(2);
 521   InitPadCOLZ();
 522   corrVtx->DrawCopy("COLZ");
 523
 524   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
 525   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
 526
 527   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
 528   canvas->Divide(2, 1);
 529
 530   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
 531   corrVtx->GetYaxis()->SetRangeUser(0, 5);
 532
 533   canvas->cd(1);
 534   InitPadCOLZ();
 535   corrTrigger->DrawCopy("COLZ");
 536
 537   canvas->cd(2);
 538   InitPadCOLZ();
 539   corrVtx->DrawCopy("COLZ");
 540
 541   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
 542   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
 543 }
 544
 545 void TriggerBias(const char* fileName = "correction_map.root")
 546 {
 547   TFile* file = TFile::Open(fileName);
 548
 549   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
 550
 551   Prepare2DPlot(corr);
 552   corr->SetTitle("Trigger bias correction");
 553
 554   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
 555   InitPadCOLZ();
 556   corr->DrawCopy("COLZ");
 557
 558   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
 559   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
 560
 561   corr->GetYaxis()->SetRangeUser(0, 5);
 562
 563   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
 564   InitPadCOLZ();
 565   corr->DrawCopy("COLZ");
 566
 567   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
 568   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
 569 }
 570
 571 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
 572 {
 573   gSystem->Load("libPWG0base");
 574
 575   TFile* file = TFile::Open(fileName);
 576   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
 577   dNdEtaCorrection->LoadHistograms(fileName, folderName);
 578
 579   TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrection()->Get1DCorrection("x");
 580   TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrection()->Get1DCorrection("y", -10, 10);
 581
 582   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
 583   canvas->Divide(2, 1);
 584
 585   canvas->cd(1);
 586   InitPad();
 587
 588   Prepare1DPlot(hist);
 589   hist->SetTitle("");
 590   hist->GetYaxis()->SetTitle("correction factor");
 591   hist->GetYaxis()->SetRangeUser(1, 1.5);
 592   hist->GetYaxis()->SetTitleOffset(1.6);
 593   hist->Draw();
 594
 595   canvas->cd(2);
 596   InitPad();
 597
 598   Prepare1DPlot(hist2);
 599   hist2->SetTitle("");
 600   hist2->GetYaxis()->SetTitle("correction factor");
 601   hist2->GetXaxis()->SetRangeUser(0, 5);
 602   hist2->GetYaxis()->SetTitleOffset(1.6);
 603   hist2->GetXaxis()->SetTitle("multiplicity");
 604   hist2->Draw();
 605
 606   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
 607   pave->SetFillColor(0);
 608   pave->AddText("|z| < 10 cm");
 609   pave->Draw();
 610
 611   canvas->SaveAs("TriggerBias1D.eps");
 612 }
 613
 614 void VtxRecon()
 615 {
 616   TFile* file = TFile::Open("correction_map.root");
 617
 618   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
 619
 620   Prepare2DPlot(corr);
 621   corr->SetTitle("Vertex reconstruction correction");
 622
 623   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
 624   InitPadCOLZ();
 625   corr->DrawCopy("COLZ");
 626
 627   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
 628   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
 629
 630   corr->GetYaxis()->SetRangeUser(0, 5);
 631
 632   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
 633   InitPadCOLZ();
 634   corr->DrawCopy("COLZ");
 635
 636   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
 637   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
 638 }
 639
 640 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
 641 {
 642   gSystem->Load("libPWG0base");
 643
 644   TFile* file = TFile::Open(fileName);
 645   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
 646   dNdEtaCorrection->LoadHistograms(fileName, folderName);
 647
 648   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->Get1DCorrection("x");
 649   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->Get1DCorrection("y", -10, 10);
 650
 651   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
 652   canvas->Divide(2, 1);
 653
 654   canvas->cd(1);
 655   InitPad();
 656
 657   Prepare1DPlot(hist);
 658   hist->SetTitle("");
 659   hist->GetYaxis()->SetTitle("correction factor");
 660   hist->GetYaxis()->SetRangeUser(1, 1.8);
 661   hist->GetYaxis()->SetTitleOffset(1.6);
 662   hist->DrawCopy();
 663
 664   canvas->cd(2);
 665   InitPad();
 666
 667   Prepare1DPlot(hist2);
 668   hist2->SetTitle("");
 669   hist2->GetYaxis()->SetTitle("correction factor");
 670   hist2->GetXaxis()->SetRangeUser(0, 20);
 671   hist2->GetYaxis()->SetTitleOffset(1.6);
 672   hist2->GetXaxis()->SetTitle("multiplicity");
 673   hist2->Draw();
 674
 675   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
 676   pave->SetFillColor(0);
 677   pave->AddText("|z| < 10 cm");
 678   pave->Draw();
 679
 680   canvas->SaveAs("VtxRecon1D.eps");
 681
 682   for (Int_t i=1; i<=hist->GetNbinsX() / 2; ++i)
 683     if (hist->GetBinContent(hist->GetNbinsX() + 1 - i) != 0)
 684       hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinContent(hist->GetNbinsX() + 1 - i));
 685
 686   new TCanvas;
 687   hist->GetXaxis()->SetRange(1, hist->GetNbinsX() / 2);
 688   hist->GetYaxis()->SetRangeUser(0.8, 1.2);
 689   hist->DrawCopy();
 690 }
 691
 692 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
 693 {
 694   gSystem->Load("libPWG0base");
 695
 696   TFile::Open(fileName);
 697   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 698   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
 699
 700   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 701   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 702
 703   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 704   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 705
 706   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 707   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
 708   gene->GetXaxis()->SetRangeUser(-10, 10);
 709   meas->GetXaxis()->SetRangeUser(-10, 10);
 710
 711   Float_t eff1 = gene->Integral() / meas->Integral();
 712   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
 713
 714   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
 715
 716   gene->GetZaxis()->SetRangeUser(0.3, 10);
 717   meas->GetZaxis()->SetRangeUser(0.3, 10);
 718
 719   Float_t eff2 = gene->Integral() / meas->Integral();
 720   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
 721
 722   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
 723
 724   gene->GetZaxis()->SetRangeUser(0.3, 1);
 725   meas->GetZaxis()->SetRangeUser(0.3, 1);
 726
 727   Float_t eff3 = gene->Integral() / meas->Integral();
 728   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
 729
 730   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
 731 }
 732
 733 void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 10)
 734 {
 735   TFile::Open(fileName);
 736   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
 737   dNdEtaCorrection->LoadHistograms(fileName, folderName);
 738
 739   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 740   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 741
 742   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
 743   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
 744   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 745   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
 746   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kTRUE);
 747   gene->GetYaxis()->SetRange(0, 0);
 748   meas->GetYaxis()->SetRange(0, 0);
 749
 750   gene->GetXaxis()->SetRangeUser(-10, 10);
 751   meas->GetXaxis()->SetRangeUser(-10, 10);
 752   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kTRUE);
 753   gene->GetZaxis()->SetRange(0, 0);
 754   meas->GetZaxis()->SetRange(0, 0);
 755
 756   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 757   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
 758   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kTRUE);
 759 }
 760
 761 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
 762 {
 763   gSystem->Load("libPWG0base");
 764
 765   Track2Particle1DCreatePlots(fileName, folder, upperPtLimit);
 766
 767   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", folder, folder)));
 768   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", folder, folder)));
 769   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", folder, folder)));
 770
 771   Prepare1DPlot(corrX);
 772   Prepare1DPlot(corrY);
 773   Prepare1DPlot(corrZ);
 774
 775   //corrX->SetTitle("a) z projection");
 776   corrY->SetTitle("a) #eta projection");
 777   corrZ->SetTitle("b) p_{T} projection");
 778
 779   corrY->GetYaxis()->SetTitle("correction factor");
 780   corrZ->GetYaxis()->SetTitle("correction factor");
 781
 782   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
 783
 784   TString canvasName;
 785   canvasName.Form("Track2Particle1D_%s", folder);
 786   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
 787   canvas->Divide(3, 1);
 788
 789   canvas->cd(1);
 790   InitPad();
 791   corrX->DrawCopy();
 792
 793   canvas->cd(2);
 794   InitPad();
 795   corrY->DrawCopy();
 796
 797   canvas->cd(3);
 798   InitPad();
 799   corrZ->DrawCopy();
 800
 801   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
 802   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
 803
 804   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
 805
 806   canvasName.Form("Track2Particle1D_%s_etapt", folder);
 807   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
 808   canvas->Divide(2, 1);
 809
 810   canvas->cd(1);
 811   InitPad();
 812   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
 813   corrY->GetYaxis()->SetRangeUser(1, 1.5);
 814   corrY->GetYaxis()->SetTitleOffset(1.5);
 815   corrY->DrawCopy();
 816   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
 817   pave->AddText("|z| < 10 cm");
 818   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
 819   pave->Draw();
 820
 821   canvas->cd(2);
 822   InitPad();
 823   gPad->SetLogx();
 824   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
 825   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
 826   corrZ->GetYaxis()->SetTitleOffset(1.5);
 827   corrZ->DrawCopy();
 828   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
 829   pave->AddText("|z| < 10 cm");
 830   pave->AddText("|#eta| < 0.8");
 831   pave->Draw();
 832
 833   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
 834   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
 835 }
 836
 837 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
 838 {
 839   gSystem->Load("libPWG0base");
 840
 841   // particle type
 842   for (Int_t particle=0; particle<4; ++particle)
 843   {
 844     TString dirName;
 845     dirName.Form("correction_%d", particle);
 846     Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
 847
 848     TString tmpx, tmpy, tmpz;
 849     tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
 850     tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
 851     tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
 852
 853     TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
 854     TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
 855     TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
 856
 857     Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
 858
 859     TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
 860     TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
 861     TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
 862
 863     posX->Divide(negX);
 864     posY->Divide(negY);
 865     posZ->Divide(negZ);
 866
 867     Prepare1DPlot(posX);
 868     Prepare1DPlot(posY);
 869     Prepare1DPlot(posZ);
 870
 871     Float_t min = 0.8;
 872     Float_t max = 1.2;
 873
 874     posX->SetMinimum(min);
 875     posX->SetMaximum(max);
 876     posY->SetMinimum(min);
 877     posY->SetMaximum(max);
 878     posZ->SetMinimum(min);
 879     posZ->SetMaximum(max);
 880
 881     posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
 882
 883     posX->GetYaxis()->SetTitleOffset(1.7);
 884     posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
 885     posY->GetYaxis()->SetTitleOffset(1.7);
 886     posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
 887     posZ->GetYaxis()->SetTitleOffset(1.7);
 888     posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
 889
 890     posZ->GetXaxis()->SetRangeUser(0, 1);
 891
 892     TString canvasName;
 893     canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
 894
 895     TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
 896     canvas->Divide(3, 1);
 897
 898     canvas->cd(1);
 899     InitPad();
 900     posX->DrawCopy();
 901
 902     canvas->cd(2);
 903     InitPad();
 904     posY->DrawCopy();
 905
 906     canvas->cd(3);
 907     InitPad();
 908     posZ->DrawCopy();
 909
 910     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
 911     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 912   }
 913 }
 914
 915 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
 916 {
 917   TFile::Open(fileName);
 918   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 919   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
 920
 921   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 922   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 923
 924   gene->GetZaxis()->SetRangeUser(0.3, 10);
 925   meas->GetZaxis()->SetRangeUser(0.3, 10);
 926   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
 927   gene->GetZaxis()->SetRange(0, 0);
 928   meas->GetZaxis()->SetRange(0, 0);
 929
 930   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 931   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
 932   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
 933   gene->GetYaxis()->SetRange(0, 0);
 934   meas->GetYaxis()->SetRange(0, 0);
 935
 936   gene->GetXaxis()->SetRangeUser(-10, 10);
 937   meas->GetXaxis()->SetRangeUser(-10, 10);
 938   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
 939   gene->GetXaxis()->SetRange(0, 0);
 940   meas->GetXaxis()->SetRange(0, 0);
 941 }
 942
 943 void Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
 944 {
 945   gSystem->Load("libPWG0base");
 946
 947   Track2Particle2DCreatePlots(fileName);
 948
 949   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_yx_div_meas_%s_nTrackToNPart_yx", folder, folder)));
 950   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_zx_div_meas_%s_nTrackToNPart_zx", folder, folder)));
 951   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_zy_div_meas_%s_nTrackToNPart_zy", folder, folder)));
 952
 953   /* this reads them from the file
 954   TH2* corrYX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx"));
 955   TH2* corrZX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx"));
 956   TH2* corrZY = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy"));*/
 957
 958   Prepare2DPlot(corrYX);
 959   Prepare2DPlot(corrZX);
 960   Prepare2DPlot(corrZY);
 961
 962   const char* title = "";
 963   corrYX->SetTitle(title);
 964   corrZX->SetTitle(title);
 965   corrZY->SetTitle(title);
 966
 967   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
 968   canvas->Divide(3, 1);
 969
 970   canvas->cd(1);
 971   InitPadCOLZ();
 972   corrYX->Draw("COLZ");
 973
 974   canvas->cd(2);
 975   InitPadCOLZ();
 976   corrZX->Draw("COLZ");
 977
 978   canvas->cd(3);
 979   InitPadCOLZ();
 980   corrZY->Draw("COLZ");
 981
 982   canvas->SaveAs(Form("Track2Particle2D_%d.gif", gMax));
 983   canvas->SaveAs(Form("Track2Particle2D_%d.eps", gMax));
 984 }
 985
 986 void CompareTrack2Particle2D()
 987 {
 988   gSystem->Load("libPWG0base");
 989
 990   Track2Particle2DCreatePlots("correction_maponly-positive.root");
 991
 992   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("pos_yx"));
 993   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("pos_zx"));
 994   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("pos_zy"));
 995
 996   Track2Particle2DCreatePlots("correction_maponly-negative.root");
 997
 998   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("neg_yx"));
 999   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("neg_zx"));
1000   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("neg_zy"));
1001
1002   posYX->Divide(negYX);
1003   posZX->Divide(negZX);
1004   posZY->Divide(negZY);
1005
1006   Prepare2DPlot(posYX);
1007   Prepare2DPlot(posZX);
1008   Prepare2DPlot(posZY);
1009
1010   Float_t min = 0.8;
1011   Float_t max = 1.2;
1012
1013   posYX->SetMinimum(min);
1014   posYX->SetMaximum(max);
1015   posZX->SetMinimum(min);
1016   posZX->SetMaximum(max);
1017   posZY->SetMinimum(min);
1018   posZY->SetMaximum(max);
1019
1020   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1021   canvas->Divide(3, 1);
1022
1023   canvas->cd(1);
1024   InitPadCOLZ();
1025   posYX->Draw("COLZ");
1026
1027   canvas->cd(2);
1028   InitPadCOLZ();
1029   posZX->Draw("COLZ");
1030
1031   canvas->cd(3);
1032   InitPadCOLZ();
1033   posZY->Draw("COLZ");
1034
1035   canvas->SaveAs("CompareTrack2Particle2D.gif");
1036   canvas->SaveAs("CompareTrack2Particle2D.eps");
1037 }
1038
1039 void Track2Particle3D()
1040 {
1041   // get left margin proper
1042
1043   TFile* file = TFile::Open("correction_map.root");
1044
1045   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1046
1047   corr->SetTitle("Correction Factor");
1048   SetRanges(corr->GetZaxis());
1049
1050   Prepare3DPlot(corr);
1051
1052   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1053   canvas->SetTheta(29.428);
1054   canvas->SetPhi(16.5726);
1055
1056   corr->Draw();
1057
1058   canvas->SaveAs("Track2Particle3D.gif");
1059   canvas->SaveAs("Track2Particle3D.eps");
1060 }
1061
1062 void Track2Particle3DAll()
1063 {
1064   TFile* file = TFile::Open("correction_map.root");
1065
1066   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1067   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1068   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1069
1070   gene->SetTitle("Generated Particles");
1071   meas->SetTitle("Measured Tracks");
1072   corr->SetTitle("Correction Factor");
1073
1074   Prepare3DPlot(gene);
1075   Prepare3DPlot(meas);
1076   Prepare3DPlot(corr);
1077
1078   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1079   canvas->Divide(3, 1);
1080
1081   canvas->cd(1);
1082   InitPad();
1083   gene->Draw();
1084
1085   canvas->cd(2);
1086   meas->Draw();
1087
1088   canvas->cd(3);
1089   corr->Draw();
1090
1091   canvas->SaveAs("Track2Particle3DAll.gif");
1092   canvas->SaveAs("Track2Particle3DAll.eps");
1093 }
1094
1095 void MultiplicityMC(Int_t xRangeMax = 50)
1096 {
1097   TFile* file = TFile::Open("multiplicityMC.root");
1098
1099   if (!file)
1100   {
1101     printf("multiplicityMC.root could not be opened.\n");
1102     return;
1103   }
1104
1105   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1106   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1107   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1108
1109   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1110   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1111   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1112   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1113   {
1114     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1115     proj->Fit("gaus", "0");
1116     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1117     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1118
1119     continue;
1120
1121     // draw for debugging
1122     new TCanvas;
1123     proj->DrawCopy();
1124     proj->GetFunction("gaus")->DrawCopy("SAME");
1125   }
1126
1127   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1128
1129   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1130   {
1131     Float_t mean = correction->GetBinContent(i);
1132     Float_t width = correctionWidth->GetBinContent(i);
1133
1134     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1135     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1136     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1137
1138     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1139     {
1140       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1141     }
1142   }
1143
1144   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1145   fMultiplicityESDCorrectedRebinned->Rebin(10);
1146   fMultiplicityESDCorrectedRebinned->Scale(0.1);
1147
1148   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1149   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1150   ratio->Divide(fMultiplicityMC);
1151
1152   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1153   ratio2->Divide(fMultiplicityMC);
1154
1155   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1156   canvas->Divide(3, 2);
1157
1158   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1159   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1160   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1161   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1162   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1163   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1164   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1165
1166   canvas->cd(1); //InitPad();
1167   fMultiplicityESD->Draw();
1168   fMultiplicityMC->SetLineColor(2);
1169   fMultiplicityMC->Draw("SAME");
1170
1171   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1172   legend->AddEntry(fMultiplicityESD, "ESD");
1173   legend->AddEntry(fMultiplicityMC, "MC");
1174   legend->Draw();
1175
1176   canvas->cd(2);
1177   fCorrelation->Draw("COLZ");
1178
1179   canvas->cd(3);
1180   correction->Draw();
1181   //correction->Fit("pol1");
1182   correctionWidth->SetLineColor(2);
1183   correctionWidth->Draw("SAME");
1184
1185   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1186   legend->AddEntry(correction, "#bar{x}");
1187   legend->AddEntry(correctionWidth, "#sigma");
1188   legend->Draw();
1189
1190   canvas->cd(4);
1191   ratio->Draw();
1192
1193   ratio2->SetLineColor(2);
1194   ratio2->Draw("SAME");
1195
1196   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1197   legend->AddEntry(ratio, "uncorrected");
1198   legend->AddEntry(ratio2, "corrected");
1199   legend->Draw();
1200
1201   canvas->cd(5);
1202   fMultiplicityESDCorrected->SetLineColor(kBlue);
1203   fMultiplicityESDCorrected->Draw();
1204   fMultiplicityMC->Draw("SAME");
1205   fMultiplicityESD->Draw("SAME");
1206
1207   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1208   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1209   legend->AddEntry(fMultiplicityMC, "MC");
1210   legend->AddEntry(fMultiplicityESD, "ESD");
1211   legend->Draw();
1212
1213   canvas->cd(6);
1214   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
1215   fMultiplicityESDCorrectedRebinned->Draw();
1216   fMultiplicityMC->Draw("SAME");
1217
1218   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1219   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1220   legend->AddEntry(fMultiplicityMC, "MC");
1221   legend->Draw();
1222
1223   canvas->SaveAs("MultiplicityMC.gif");
1224 }
1225
1226 void MultiplicityESD()
1227 {
1228   TFile* file = TFile::Open("multiplicityESD.root");
1229
1230   if (!file)
1231   {
1232     printf("multiplicityESD.root could not be opened.\n");
1233     return;
1234   }
1235
1236   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1237
1238   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1239
1240   fMultiplicityESD->Draw();
1241 }
1242
1243 void drawPlots(Int_t max)
1244 {
1245   gMax = max;
1246
1247   ptCutoff();
1248   TriggerBias();
1249   VtxRecon();
1250   Track2Particle2D();
1251   Track2Particle3D();
1252   ptSpectrum();
1253   dNdEta();
1254 }
1255
1256 void drawPlots()
1257 {
1258   drawPlots(5);
1259   drawPlots(2);
1260 }