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 loadlibs()
  24 {
  25   gSystem->Load("libANALYSIS");
  26   gSystem->Load("libPWG0base");
  27 }
  28
  29 void SetRanges(TAxis* axis)
  30 {
  31   if (strcmp(axis->GetTitle(), "#eta") == 0)
  32     axis->SetRangeUser(-1.7999, 1.7999);
  33   if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
  34     axis->SetRangeUser(0, 9.9999);
  35   if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
  36     axis->SetRangeUser(-15, 14.9999);
  37   if (strcmp(axis->GetTitle(), "Ntracks") == 0)
  38     axis->SetRangeUser(0, 99.9999);
  39 }
  40
  41 void SetRanges(TH1* hist)
  42 {
  43   SetRanges(hist->GetXaxis());
  44   SetRanges(hist->GetYaxis());
  45   SetRanges(hist->GetZaxis());
  46 }
  47
  48
  49 void Prepare3DPlot(TH3* hist)
  50 {
  51   hist->GetXaxis()->SetTitleOffset(1.5);
  52   hist->GetYaxis()->SetTitleOffset(1.5);
  53   hist->GetZaxis()->SetTitleOffset(1.5);
  54
  55   hist->SetStats(kFALSE);
  56 }
  57
  58 void Prepare2DPlot(TH2* hist)
  59 {
  60   hist->SetStats(kFALSE);
  61   hist->GetYaxis()->SetTitleOffset(1.4);
  62
  63   hist->SetMinimum(0);
  64   hist->SetMaximum(gMax);
  65
  66   SetRanges(hist);
  67 }
  68
  69 void Prepare1DPlot(TH1* hist)
  70 {
  71   hist->SetLineWidth(2);
  72   hist->SetStats(kFALSE);
  73
  74   hist->GetXaxis()->SetLabelOffset(0.02);
  75   hist->GetXaxis()->SetTitleOffset(1.3);
  76   hist->GetYaxis()->SetTitleOffset(1.3);
  77
  78   SetRanges(hist);
  79 }
  80
  81 void InitPad()
  82 {
  83   gPad->Range(0, 0, 1, 1);
  84   gPad->SetLeftMargin(0.15);
  85   //gPad->SetRightMargin(0.05);
  86   //gPad->SetTopMargin(0.13);
  87   gPad->SetBottomMargin(0.12);
  88
  89   gPad->SetGridx();
  90   gPad->SetGridy();
  91 }
  92
  93 void InitPadCOLZ()
  94 {
  95   gPad->Range(0, 0, 1, 1);
  96   gPad->SetRightMargin(0.15);
  97   gPad->SetLeftMargin(0.12);
  98
  99   gPad->SetGridx();
 100   gPad->SetGridy();
 101 }
 102
 103 // --- end of helpers --- begin functions ---
 104
 105 void DrawOverview(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
 106 {
 107   loadlibs();
 108   if (!TFile::Open(fileName))
 109     return;
 110
 111   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
 112   if (!dNdEtaCorrection->LoadHistograms())
 113     return;
 114
 115   dNdEtaCorrection->Finish();
 116
 117   dNdEtaCorrection->DrawOverview();
 118 }
 119
 120 void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
 121 {
 122   gSystem->Load("libPWG0base");
 123
 124   TFile::Open(esdFile);
 125   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 126   fdNdEtaAnalysisESD->LoadHistograms();
 127
 128   TFile::Open(mcFile);
 129   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
 130   fdNdEtaAnalysisMC->LoadHistograms();
 131   //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 132
 133   for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
 134     fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
 135
 136   fdNdEtaAnalysisESD->DrawHistograms();
 137 }
 138
 139 void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 140 {
 141   gSystem->Load("libPWG0base");
 142
 143   const char* ESDfolder = 0;
 144
 145   if (plot == 0) // all
 146     ESDfolder = "dndeta";
 147   else if (plot == 1) // mb
 148     ESDfolder = "dndeta_mb";
 149   else if (plot == 2) // mb vtx
 150     ESDfolder = "dndeta_mbvtx";
 151
 152   TFile::Open("analysis_esd.root");
 153   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 154   fdNdEtaAnalysisESD->LoadHistograms();
 155
 156   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 157   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 158
 159   TH2F* hist = 0;
 160
 161   if (plot == 0) // all
 162     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
 163   else if (plot == 1) // mb
 164     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
 165   else if (plot == 2) // mb vtx
 166     hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
 167
 168   TH1* proj = hist->ProjectionX();
 169
 170   TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
 171   for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
 172   {
 173     Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
 174     if (value != 0)
 175     {
 176       printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
 177       vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
 178     }
 179   }
 180
 181   new TCanvas;
 182   vertex->DrawCopy();
 183 }
 184
 185 void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 186 {
 187   gSystem->Load("libPWG0base");
 188
 189   TFile::Open("analysis_esd.root");
 190   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
 191   fdNdEtaAnalysisESD->LoadHistograms();
 192
 193   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 194   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 195
 196   //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
 197   //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 198
 199   TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
 200   TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 201
 202   TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
 203
 204   new TCanvas;
 205   histESD->Draw();
 206
 207   new TCanvas;
 208   histCorr->Draw();
 209
 210   for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
 211     for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
 212       for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
 213       {
 214         Float_t value1 = histESD->GetBinContent(x, y, z);
 215         Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
 216
 217         if (value2 > 0 && value1 > 0)
 218         {
 219           printf("%f %f %f\n", value1, value2, value1 / value2);
 220           diff->Fill(value1 / value2);
 221         }
 222       }
 223
 224   new TCanvas;
 225   diff->Draw();
 226 }
 227
 228 Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
 229 {
 230   Double_t avgMC = 0;
 231   Double_t avgESD = 0;
 232   for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
 233   {
 234     avgMC += histMC->GetBinContent(bin);
 235     avgESD += histESD->GetBinContent(bin);
 236   }
 237   Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
 238
 239   avgMC /= nBins;
 240   avgESD /= nBins;
 241
 242   // deviation when integrate in |eta| < 0.8 between mc and esd
 243   Double_t diffFullRange = (avgMC - avgESD) / avgMC;
 244
 245   Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
 246
 247   return diffFullRange;
 248 }
 249
 250 void dNdEtaNoResolution()
 251 {
 252   loadlibs();
 253
 254   TFile::Open("correction_map.root");
 255
 256   const char* correctionMapFolder = "dndeta_correction";
 257   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 258   dNdEtaCorrection->LoadHistograms();
 259   dNdEtaCorrection->GetTrack2ParticleCorrection()->PrintInfo(0.3);
 260
 261   TFile::Open("analysis_mc.root");
 262   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
 263   fdNdEtaAnalysis->LoadHistograms();
 264   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, "ESD (no resolution effect) -> MB with vertex");
 265   fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0)->SetMarkerStyle(21);
 266
 267   TFile::Open("analysis_mc.root");
 268   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
 269   fdNdEtaAnalysisMC->LoadHistograms();
 270   fdNdEtaAnalysisMC->Finish(0, 0, AlidNdEtaCorrection::kNone, "MC: MB with vertex");
 271
 272   DrawdNdEtaRatio(fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0), fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(0), "MB with vertex (no resolution effect)", 3);
 273 }
 274
 275 TH1* GetMCHist(const char* folder, Float_t ptCut, const char* tag)
 276 {
 277   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(folder, folder);
 278   fdNdEtaAnalysis->LoadHistograms();
 279   fdNdEtaAnalysis->Finish(0, ptCut, AlidNdEtaCorrection::kNone, tag);
 280   return fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 281 }
 282
 283 void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 284 {
 285   TFile* file = TFile::Open("analysis_esd.root");
 286   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
 287   TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
 288   TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
 289   TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
 290   TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
 291   TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
 292   TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
 293   TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
 294
 295   Prepare1DPlot(histESD);
 296   Prepare1DPlot(histESDnsd);
 297   Prepare1DPlot(histESDMB);
 298   Prepare1DPlot(histESDMBVtx);
 299
 300   Prepare1DPlot(histESDNoPt);
 301   Prepare1DPlot(histESDMBNoPt);
 302   Prepare1DPlot(histESDMBVtxNoPt);
 303   Prepare1DPlot(histESDMBTracksNoPt);
 304
 305   histESD->SetLineWidth(0);
 306   histESDnsd->SetLineWidth(0);
 307   histESDMB->SetLineWidth(0);
 308   histESDMBVtx->SetLineWidth(0);
 309
 310   histESDNoPt->SetLineWidth(0);
 311   histESDMBNoPt->SetLineWidth(0);
 312   histESDMBVtxNoPt->SetLineWidth(0);
 313
 314   histESD->SetMarkerColor(1);
 315   histESDnsd->SetMarkerColor(6);
 316   histESDMB->SetMarkerColor(2);
 317   histESDMBVtx->SetMarkerColor(3);
 318
 319   histESD->SetLineColor(1);
 320   histESDnsd->SetLineColor(6);
 321   histESDMB->SetLineColor(2);
 322   histESDMBVtx->SetLineColor(3);
 323
 324   histESDNoPt->SetMarkerColor(1);
 325   histESDMBNoPt->SetMarkerColor(2);
 326   histESDMBVtxNoPt->SetMarkerColor(3);
 327   histESDMBTracksNoPt->SetMarkerColor(4);
 328
 329   histESD->SetMarkerStyle(20);
 330   histESDnsd->SetMarkerStyle(29);
 331   histESDMB->SetMarkerStyle(21);
 332   histESDMBVtx->SetMarkerStyle(22);
 333
 334   histESDNoPt->SetMarkerStyle(20);
 335   histESDMBNoPt->SetMarkerStyle(21);
 336   histESDMBVtxNoPt->SetMarkerStyle(22);
 337   histESDMBTracksNoPt->SetMarkerStyle(23);
 338
 339   //Float_t etaLimit = 1.2999;
 340   Float_t etaLimit = 2.41;
 341   Float_t etaPlotLimit = 2.6;
 342
 343   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 344   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 345   histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 346   histESDnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 347
 348   histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 349   histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 350   histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 351   histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 352
 353   Float_t max = TMath::Max(histESDMBVtx->GetMaximum(), histESDMB->GetMaximum());
 354   max = TMath::Max(max, histESD->GetMaximum());
 355
 356   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, max * 1.1);
 357   Prepare1DPlot(dummy);
 358   dummy->SetStats(kFALSE);
 359   dummy->SetXTitle("#eta");
 360   dummy->SetYTitle("dN_{ch}/d#eta");
 361   dummy->GetYaxis()->SetTitleOffset(1);
 362
 363   TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
 364
 365   dummy->DrawCopy();
 366   histESDMBVtx->Draw("SAME");
 367   histESDMB->Draw("SAME");
 368   histESD->Draw("SAME");
 369
 370   if (save)
 371   {
 372     canvas->SaveAs("dNdEta1.gif");
 373     canvas->SaveAs("dNdEta1.eps");
 374   }
 375
 376   if (onlyESD)
 377     return;
 378
 379   loadlibs();
 380
 381   TFile* file2 = TFile::Open("analysis_mc.root");
 382
 383   TH1* histMC =            (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
 384   TH1* histMCTr =          (TH1*) GetMCHist("dndetaTr", -1, "MC: minimum bias")->Clone("histMCTr");
 385   TH1* histMCTrVtx =       (TH1*) GetMCHist("dndetaTrVtx", -1, "MC: MB with trigger")->Clone("histMCTrVtx");
 386   TH1* histMCnsd =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
 387
 388   TH1* histMCPtCut =       (TH1*) GetMCHist("dndeta", 0.3, "MC: full inelastic, pt cut")->Clone("histMCPtCut");
 389   TH1* histMCTrPtCut =     (TH1*) GetMCHist("dndetaTr", 0.3, "MC: minimum bias, pt cut")->Clone("histMCTrPtCut");
 390   TH1* histMCTrVtxPtCut =  (TH1*) GetMCHist("dndetaTrVtx", 0.3, "MC: MB with trigger, pt cut")->Clone("histMCTrVtxPtCut");
 391   TH1* histMCTracksPtCut = (TH1*) GetMCHist("dndetaTracks", 0.3, "MC: Tracks w/o resolution effect, pt cut")->Clone("histMCTracksPtCut");
 392
 393   Prepare1DPlot(histMC);
 394   Prepare1DPlot(histMCnsd);
 395   Prepare1DPlot(histMCTr);
 396   Prepare1DPlot(histMCTrVtx);
 397
 398   Prepare1DPlot(histMCPtCut);
 399   Prepare1DPlot(histMCTrPtCut);
 400   Prepare1DPlot(histMCTrVtxPtCut);
 401   Prepare1DPlot(histMCTracksPtCut);
 402
 403   histMC->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 404   histMCnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 405   histMCTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 406   histMCTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 407
 408   histMCPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 409   histMCTrPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 410   histMCTrVtxPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 411   histMCTracksPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 412
 413   histMC->SetLineColor(1);
 414   histMCnsd->SetLineColor(6);
 415   histMCTr->SetLineColor(2);
 416   histMCTrVtx->SetLineColor(3);
 417
 418   histMCPtCut->SetLineColor(1);
 419   histMCTrPtCut->SetLineColor(2);
 420   histMCTrVtxPtCut->SetLineColor(3);
 421   if (histMCTracksPtCut)
 422     histMCTracksPtCut->SetLineColor(4);
 423
 424   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
 425
 426   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
 427   dummy2->GetYaxis()->SetRangeUser(0, max * 1.1);
 428
 429   dummy2->DrawCopy();
 430   histMC->Draw("SAME");
 431   histMCnsd->Draw("SAME");
 432   histMCTr->Draw("SAME");
 433   histMCTrVtx->Draw("SAME");
 434   histESD->Draw("SAME");
 435   histESDnsd->Draw("SAME");
 436   histESDMB->Draw("SAME");
 437   histESDMBVtx->Draw("SAME");
 438   histESDNoPt->Draw("SAME");
 439   histESDMBNoPt->Draw("SAME");
 440   histESDMBVtxNoPt->Draw("SAME");
 441   histESDMBTracksNoPt->Draw("SAME");
 442   histMCPtCut->Draw("SAME");
 443   histMCTrPtCut->Draw("SAME");
 444   histMCTrVtxPtCut->Draw("SAME");
 445   if (histMCTracksPtCut)
 446     histMCTracksPtCut->Draw("SAME");
 447
 448   if (save)
 449   {
 450     canvas2->SaveAs("dNdEta2.gif");
 451     canvas2->SaveAs("dNdEta2.eps");
 452   }
 453
 454   DrawdNdEtaRatio(histESD, histMC, "full_inelastic", etaPlotLimit);
 455   DrawdNdEtaRatio(histESDMB, histMCTr, "triggered", etaPlotLimit);
 456   DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit);
 457
 458   new TCanvas;
 459   dummy2->DrawCopy();
 460   histMCnsd->Draw("SAME");
 461   histESDnsd->Draw("SAME");
 462
 463   TH1* ratio = (TH1*) histMC->Clone("ratio");
 464   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
 465
 466   ratio->Divide(histESD);
 467   ratioNoPt->Divide(histESDNoPt);
 468
 469   ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 470
 471   ratio->SetLineColor(1);
 472   ratioNoPt->SetLineColor(2);
 473
 474   Double_t average = 0;       // average deviation from 1 in ratio (depends on the number of bins if statistical)
 475   for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
 476     average += TMath::Abs(ratio->GetBinContent(bin) - 1);
 477   Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
 478   average /= nBins;
 479   Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
 480
 481   PrintIntegratedDeviation(histMC, histESD, "all events");
 482   PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
 483   PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
 484   PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
 485   PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
 486   PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
 487
 488   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 600);
 489   canvas3->Range(0, 0, 1, 1);
 490   //canvas3->Divide(1, 2, 0, 0);
 491
 492   //canvas3->cd(1);
 493   TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
 494   pad1->Draw();
 495
 496   TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
 497   pad2->Draw();
 498
 499   pad1->SetRightMargin(0.05);
 500   pad2->SetRightMargin(0.05);
 501
 502   // no border between them
 503   pad1->SetBottomMargin(0);
 504   pad2->SetTopMargin(0);
 505
 506   pad1->cd();
 507
 508   TLegend* legend = new TLegend(0.4, 0.05, 0.65, 0.3);
 509   legend->SetFillColor(0);
 510   legend->AddEntry(histESDMBVtx, "triggered, vertex");
 511   legend->AddEntry(histESDMB, "triggered");
 512   legend->AddEntry(histESD, "all events");
 513   legend->AddEntry(histMC, "MC prediction");
 514
 515   dummy->GetXaxis()->SetLabelSize(0.06);
 516   dummy->GetYaxis()->SetLabelSize(0.06);
 517   dummy->GetXaxis()->SetTitleSize(0.06);
 518   dummy->GetYaxis()->SetTitleSize(0.06);
 519   dummy->GetYaxis()->SetTitleOffset(0.7);
 520   dummy->DrawCopy();
 521   histESDMBVtx->Draw("SAME");
 522   histESDMB->Draw("SAME");
 523   histESD->Draw("SAME");
 524   histMC->Draw("SAME");
 525
 526   legend->Draw();
 527
 528   pad2->cd();
 529   pad2->SetBottomMargin(0.15);
 530
 531   Float_t minR = 0.9; //TMath::Min(0.961, ratio->GetMinimum() * 0.95);
 532   Float_t maxR = 1.1; //TMath::Max(1.049, ratio->GetMaximum() * 1.05);
 533
 534   TH1F dummy3("dummy3", ";#eta;Ratio: MC / ESD", 100, -etaPlotLimit, etaPlotLimit);
 535   dummy3.SetStats(kFALSE);
 536   for (Int_t i=1; i<=100; ++i)
 537     dummy3.SetBinContent(i, 1);
 538   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
 539   dummy3.SetLineWidth(2);
 540   dummy3.GetXaxis()->SetLabelSize(0.06);
 541   dummy3.GetYaxis()->SetLabelSize(0.06);
 542   dummy3.GetXaxis()->SetTitleSize(0.06);
 543   dummy3.GetYaxis()->SetTitleSize(0.06);
 544   dummy3.GetYaxis()->SetTitleOffset(0.7);
 545   dummy3.DrawCopy();
 546
 547   ratio->Draw("SAME");
 548
 549   //pad2->Draw();
 550
 551   canvas3->Modified();
 552
 553   if (save)
 554   {
 555     canvas3->SaveAs("dNdEta.gif");
 556     canvas3->SaveAs("dNdEta.eps");
 557   }
 558
 559   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
 560
 561   ratio->Draw();
 562   ratioNoPt->Draw("SAME");
 563
 564   TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
 565   legend->SetFillColor(0);
 566   legend->AddEntry(ratio, "mc/esd");
 567   legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
 568   legend->Draw();
 569 }
 570
 571 void DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
 572 {
 573   TCanvas* canvas3 = new TCanvas(name, name, 700, 600);
 574   canvas3->Range(0, 0, 1, 1);
 575
 576   TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
 577   pad1->Draw();
 578
 579   TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
 580   pad2->Draw();
 581
 582   pad1->SetRightMargin(0.05);
 583   pad2->SetRightMargin(0.05);
 584
 585   // no border between them
 586   pad1->SetBottomMargin(0);
 587   pad2->SetTopMargin(0);
 588
 589   pad1->cd();
 590
 591   TLegend* legend = new TLegend(0.4, 0.05, 0.65, 0.3);
 592   legend->SetFillColor(0);
 593   legend->AddEntry(corr, "corrected");
 594   legend->AddEntry(mc, "MC prediction");
 595
 596   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, corr->GetMaximum() * 1.1);
 597   Prepare1DPlot(dummy);
 598   dummy->SetStats(kFALSE);
 599   dummy->SetXTitle("#eta");
 600   dummy->SetYTitle("dN_{ch}/d#eta");
 601   dummy->GetYaxis()->SetTitleOffset(1);
 602
 603   dummy->GetXaxis()->SetLabelSize(0.06);
 604   dummy->GetYaxis()->SetLabelSize(0.06);
 605   dummy->GetXaxis()->SetTitleSize(0.06);
 606   dummy->GetYaxis()->SetTitleSize(0.06);
 607   dummy->GetYaxis()->SetTitleOffset(0.7);
 608   dummy->DrawCopy();
 609
 610   corr->Draw("SAME");
 611   mc->Draw("SAME");
 612
 613   legend->Draw();
 614
 615   pad2->cd();
 616   pad2->SetBottomMargin(0.15);
 617
 618   TH1* ratio = (TH1*) mc->Clone("ratio");
 619   ratio->Divide(corr);
 620
 621   Float_t minR = TMath::Min(0.961, ratio->GetMinimum() * 0.95);
 622   Float_t maxR = TMath::Max(1.049, ratio->GetMaximum() * 1.05);
 623
 624   TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
 625   dummy3.SetStats(kFALSE);
 626   for (Int_t i=1; i<=100; ++i)
 627         dummy3.SetBinContent(i, 1);
 628   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
 629   dummy3.SetLineWidth(2);
 630   dummy3.GetXaxis()->SetLabelSize(0.06);
 631   dummy3.GetYaxis()->SetLabelSize(0.06);
 632   dummy3.GetXaxis()->SetTitleSize(0.06);
 633   dummy3.GetYaxis()->SetTitleSize(0.06);
 634   dummy3.GetYaxis()->SetTitleOffset(0.7);
 635   dummy3.DrawCopy();
 636
 637   ratio->Draw("SAME");
 638
 639   canvas3->Modified();
 640 }
 641
 642 void ptSpectrum()
 643 {
 644   TFile* file = TFile::Open("analysis_esd.root");
 645   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
 646
 647   TFile* file2 = TFile::Open("analysis_mc.root");
 648   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
 649
 650   TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
 651   InitPad();
 652   gPad->SetLogy();
 653
 654   Prepare1DPlot(histMC);
 655   Prepare1DPlot(histESD);
 656
 657   histESD->SetTitle("");
 658   histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 659   histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
 660
 661   histMC->SetLineColor(kBlue);
 662   histESD->SetLineColor(kRed);
 663
 664   histESD->GetYaxis()->SetTitleOffset(1.5);
 665   histESD->GetXaxis()->SetRangeUser(0, 4.9999);
 666
 667   histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
 668
 669   histESD->Draw();
 670   histMC->Draw("SAME");
 671
 672   canvas->SaveAs("ptSpectrum.gif");
 673   canvas->SaveAs("ptSpectrum.eps");
 674 }
 675
 676 void ptCutoff()
 677 {
 678   gSystem->Load("libPWG0base");
 679
 680   TFile::Open("correction_map.root");
 681   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 682   dNdEtaCorrection->LoadHistograms();
 683
 684   dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, -100, kTRUE);
 685
 686   TH1* hist = dynamic_cast<TH1*> (gROOT->FindObject("generated_pt")->Clone("ptcutoff"));
 687
 688   hist->GetXaxis()->SetRangeUser(0, 0.9999);
 689   hist->SetMinimum(0);
 690
 691   hist->SetTitle("Generated Particles");
 692   Prepare1DPlot(hist);
 693
 694   TCanvas* canvas = new TCanvas("ptCutoff", "ptCutoff", 700, 500);
 695   hist->DrawCopy();
 696
 697   TLine* line = new TLine(0.3, 0 - hist->GetMaximum() * 0, 0.3, hist->GetMaximum() * 1.1);
 698   line->SetLineWidth(3);
 699   line->SetLineColor(kRed);
 700   line->Draw();
 701
 702   canvas->SaveAs("ptCutoff.gif");
 703   canvas->SaveAs("ptCutoff.eps");
 704
 705   TH1F* factor = new TH1F("factor", ";#eta;correction factor", 20, -1, 1.000001);
 706   factor->SetLineWidth(2);
 707   for (Float_t eta = -0.95; eta<1; eta += 0.1)
 708     factor->Fill(eta, 1.0 / dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, eta, kFALSE));
 709
 710   TCanvas* canvas = new TCanvas("ptCutoff_factor", "ptCutoff_factor", 700, 500);
 711   InitPad();
 712
 713   Prepare1DPlot(factor);
 714   factor->GetYaxis()->SetRangeUser(1, 2);
 715   factor->GetYaxis()->SetTitleOffset(1);
 716   factor->Draw();
 717
 718   canvas->SaveAs("ptCutoff_factor.eps");
 719 }
 720
 721 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
 722 {
 723   gSystem->Load("libPWG0base");
 724
 725   TFile::Open(fileName);
 726   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 727   dNdEtaCorrection->LoadHistograms();
 728
 729   TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
 730   TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
 731
 732   Prepare2DPlot(corrTrigger);
 733   corrTrigger->SetTitle("b) Trigger bias correction");
 734
 735   Prepare2DPlot(corrVtx);
 736   corrVtx->SetTitle("a) Vertex reconstruction correction");
 737
 738   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
 739   corrVtx->GetYaxis()->SetTitle("Multiplicity");
 740
 741   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
 742   canvas->Divide(2, 1);
 743
 744   canvas->cd(1);
 745   InitPadCOLZ();
 746   corrVtx->DrawCopy("COLZ");
 747
 748   canvas->cd(2);
 749   InitPadCOLZ();
 750   corrTrigger->DrawCopy("COLZ");
 751
 752   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
 753   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
 754
 755   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
 756   canvas->Divide(2, 1);
 757
 758   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
 759   corrVtx->GetYaxis()->SetRangeUser(0, 5);
 760
 761   canvas->cd(1);
 762   InitPadCOLZ();
 763   corrVtx->DrawCopy("COLZ");
 764
 765   canvas->cd(2);
 766   InitPadCOLZ();
 767   corrTrigger->DrawCopy("COLZ");
 768
 769   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
 770   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
 771 }
 772
 773 void TriggerBias(const char* fileName = "correction_map.root")
 774 {
 775   TFile* file = TFile::Open(fileName);
 776
 777   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
 778
 779   Prepare2DPlot(corr);
 780   corr->SetTitle("Trigger bias correction");
 781
 782   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
 783   InitPadCOLZ();
 784   corr->DrawCopy("COLZ");
 785
 786   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
 787   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
 788
 789   corr->GetYaxis()->SetRangeUser(0, 5);
 790
 791   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
 792   InitPadCOLZ();
 793   corr->DrawCopy("COLZ");
 794
 795   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
 796   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
 797 }
 798
 799 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
 800 {
 801   gSystem->Load("libPWG0base");
 802
 803   TFile* file = TFile::Open(fileName);
 804   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
 805   dNdEtaCorrection->LoadHistograms();
 806
 807   TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
 808   TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
 809
 810   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
 811   canvas->Divide(2, 1);
 812
 813   canvas->cd(1);
 814   InitPad();
 815
 816   Prepare1DPlot(hist);
 817   hist->SetTitle("");
 818   hist->GetYaxis()->SetTitle("correction factor");
 819   hist->GetYaxis()->SetRangeUser(1, 1.5);
 820   hist->GetYaxis()->SetTitleOffset(1.6);
 821   hist->Draw();
 822
 823   canvas->cd(2);
 824   InitPad();
 825
 826   Prepare1DPlot(hist2);
 827   hist2->SetTitle("");
 828   hist2->GetYaxis()->SetTitle("correction factor");
 829   hist2->GetXaxis()->SetRangeUser(0, 5);
 830   hist2->GetYaxis()->SetTitleOffset(1.6);
 831   hist2->GetXaxis()->SetTitle("multiplicity");
 832   hist2->Draw();
 833
 834   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
 835   pave->SetFillColor(0);
 836   pave->AddText("|z| < 10 cm");
 837   pave->Draw();
 838
 839   canvas->SaveAs("TriggerBias1D.eps");
 840 }
 841
 842 void VtxRecon()
 843 {
 844   TFile* file = TFile::Open("correction_map.root");
 845
 846   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
 847
 848   Prepare2DPlot(corr);
 849   corr->SetTitle("Vertex reconstruction correction");
 850
 851   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
 852   InitPadCOLZ();
 853   corr->DrawCopy("COLZ");
 854
 855   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
 856   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
 857
 858   corr->GetYaxis()->SetRangeUser(0, 5);
 859
 860   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
 861   InitPadCOLZ();
 862   corr->DrawCopy("COLZ");
 863
 864   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
 865   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
 866 }
 867
 868 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
 869 {
 870   gSystem->Load("libPWG0base");
 871
 872   TFile* file = TFile::Open(fileName);
 873   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
 874   dNdEtaCorrection->LoadHistograms();
 875
 876   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
 877   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
 878
 879   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
 880   canvas->Divide(2, 1);
 881
 882   canvas->cd(1);
 883   InitPad();
 884
 885   Prepare1DPlot(hist);
 886   hist->SetTitle("");
 887   hist->GetYaxis()->SetTitle("correction factor");
 888   hist->GetYaxis()->SetRangeUser(1, 1.8);
 889   hist->GetYaxis()->SetTitleOffset(1.6);
 890   hist->DrawCopy();
 891
 892   canvas->cd(2);
 893   InitPad();
 894
 895   Prepare1DPlot(hist2);
 896   hist2->SetTitle("");
 897   hist2->GetYaxis()->SetTitle("correction factor");
 898   hist2->GetXaxis()->SetRangeUser(0, 20);
 899   hist2->GetYaxis()->SetTitleOffset(1.6);
 900   hist2->GetXaxis()->SetTitle("multiplicity");
 901   hist2->Draw();
 902
 903   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
 904   pave->SetFillColor(0);
 905   pave->AddText("|z| < 10 cm");
 906   pave->Draw();
 907
 908   canvas->SaveAs("VtxRecon1D.eps");
 909
 910   Correction1DCreatePlots(fileName, folderName, 9.9, 2);
 911
 912   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
 913   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
 914
 915   Prepare1DPlot(corrX);
 916   Prepare1DPlot(corrZ);
 917
 918   corrX->GetYaxis()->SetTitleOffset(1.5);
 919   corrZ->GetYaxis()->SetTitleOffset(1.5);
 920
 921   corrX->SetTitle("a) z projection");
 922   corrZ->SetTitle("b) p_{T} projection");
 923
 924   corrX->GetYaxis()->SetTitle("correction factor");
 925   corrZ->GetYaxis()->SetTitle("correction factor");
 926
 927   corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
 928
 929   TString canvasName;
 930   canvasName.Form("VtxRecon1D_Track");
 931   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
 932   canvas->Divide(2, 1);
 933
 934   canvas->cd(1);
 935   InitPad();
 936   corrX->DrawCopy();
 937
 938   canvas->cd(2);
 939   InitPad();
 940   gPad->SetLogx();
 941   corrZ->Draw();
 942
 943   canvas->SaveAs("VtxRecon1D_Track.eps");
 944   canvas->SaveAs("VtxRecon1D_Track.gif");
 945 }
 946
 947 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
 948 {
 949   gSystem->Load("libPWG0base");
 950
 951   TFile::Open(fileName);
 952   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 953   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
 954
 955   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 956   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 957
 958   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 959   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 960
 961   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 962   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
 963   gene->GetXaxis()->SetRangeUser(-10, 10);
 964   meas->GetXaxis()->SetRangeUser(-10, 10);
 965
 966   Float_t eff1 = gene->Integral() / meas->Integral();
 967   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
 968
 969   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
 970
 971   gene->GetZaxis()->SetRangeUser(0.3, 10);
 972   meas->GetZaxis()->SetRangeUser(0.3, 10);
 973
 974   Float_t eff2 = gene->Integral() / meas->Integral();
 975   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
 976
 977   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
 978
 979   gene->GetZaxis()->SetRangeUser(0.3, 1);
 980   meas->GetZaxis()->SetRangeUser(0.3, 1);
 981
 982   Float_t eff3 = gene->Integral() / meas->Integral();
 983   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
 984
 985   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
 986 }
 987
 988 void Correction1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType = 0)
 989 {
 990   TFile::Open(fileName);
 991   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
 992   dNdEtaCorrection->LoadHistograms();
 993
 994   TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
 995   TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetMeasuredHistogram();
 996
 997   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
 998   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
 999   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1000   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1001   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kTRUE);
1002   gene->GetYaxis()->SetRange(0, 0);
1003   meas->GetYaxis()->SetRange(0, 0);
1004
1005   gene->GetXaxis()->SetRangeUser(-10, 10);
1006   meas->GetXaxis()->SetRangeUser(-10, 10);
1007   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kTRUE);
1008   gene->GetZaxis()->SetRange(0, 0);
1009   meas->GetZaxis()->SetRange(0, 0);
1010
1011   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1012   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1013   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kTRUE);
1014 }
1015
1016 void Correction1D(Int_t correctionType = 0, const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1017 {
1018   gSystem->Load("libPWG0base");
1019
1020   Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType);
1021
1022   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1023   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1024   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1025
1026   Prepare1DPlot(corrX);
1027   Prepare1DPlot(corrY);
1028   Prepare1DPlot(corrZ);
1029
1030   corrX->SetTitle("a) z projection");
1031   corrY->SetTitle("b) #eta projection");
1032   corrZ->SetTitle("c) p_{T} projection");
1033
1034   corrX->GetYaxis()->SetTitle("correction factor");
1035   corrY->GetYaxis()->SetTitle("correction factor");
1036   corrZ->GetYaxis()->SetTitle("correction factor");
1037
1038   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1039
1040   TString canvasName;
1041   canvasName.Form("Correction1D_%s", folder);
1042   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1043   canvas->Divide(3, 1);
1044
1045   canvas->cd(1);
1046   InitPad();
1047   corrX->DrawCopy();
1048
1049   canvas->cd(2);
1050   InitPad();
1051   corrY->Draw();
1052
1053   canvas->cd(3);
1054   InitPad();
1055   corrZ->Draw();
1056
1057   canvas->SaveAs(Form("Correction1D_%d_%s_%f.gif", correctionType, fileName, upperPtLimit));
1058   canvas->SaveAs(Form("Correction1D_%d_%s_%f.eps", correctionType, fileName, upperPtLimit));
1059 }
1060
1061 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1062 {
1063   gSystem->Load("libPWG0base");
1064
1065   Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
1066
1067   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1068   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1069   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1070
1071   Prepare1DPlot(corrX);
1072   Prepare1DPlot(corrY);
1073   Prepare1DPlot(corrZ);
1074
1075   corrX->SetTitle("a) z projection");
1076   corrY->SetTitle("a) #eta projection");
1077   corrZ->SetTitle("b) p_{T} projection");
1078
1079   corrY->GetYaxis()->SetTitle("correction factor");
1080   corrZ->GetYaxis()->SetTitle("correction factor");
1081
1082   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1083
1084   TString canvasName;
1085   canvasName.Form("Track2Particle1D_%s", folder);
1086   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1087   canvas->Divide(3, 1);
1088
1089   canvas->cd(1);
1090   InitPad();
1091   corrX->DrawCopy();
1092
1093   canvas->cd(2);
1094   InitPad();
1095   corrY->Draw();
1096
1097   canvas->cd(3);
1098   InitPad();
1099   corrZ->Draw();
1100
1101   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
1102   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
1103
1104   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
1105
1106   canvasName.Form("Track2Particle1D_%s_etapt", folder);
1107   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
1108   canvas->Divide(2, 1);
1109
1110   canvas->cd(1);
1111   InitPad();
1112   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
1113   corrY->GetYaxis()->SetRangeUser(1, 1.5);
1114   corrY->GetYaxis()->SetTitleOffset(1.5);
1115   corrY->DrawCopy();
1116   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
1117   pave->AddText("|z| < 10 cm");
1118   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
1119   pave->Draw();
1120
1121   canvas->cd(2);
1122   InitPad();
1123   gPad->SetLogx();
1124   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
1125   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
1126   corrZ->GetYaxis()->SetTitleOffset(1.5);
1127   corrZ->DrawCopy();
1128   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
1129   pave->AddText("|z| < 10 cm");
1130   pave->AddText("|#eta| < 0.8");
1131   pave->Draw();
1132
1133   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
1134   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
1135 }
1136
1137 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1138 {
1139   gSystem->Load("libPWG0base");
1140
1141   // particle type
1142   for (Int_t particle=0; particle<4; ++particle)
1143   {
1144     TString dirName;
1145     dirName.Form("correction_%d", particle);
1146     Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
1147
1148     TString tmpx, tmpy, tmpz;
1149     tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1150     tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1151     tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
1152
1153     TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1154     TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1155     TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
1156
1157     Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
1158
1159     TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1160     TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1161     TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
1162
1163     posX->Divide(negX);
1164     posY->Divide(negY);
1165     posZ->Divide(negZ);
1166
1167     Prepare1DPlot(posX);
1168     Prepare1DPlot(posY);
1169     Prepare1DPlot(posZ);
1170
1171     Float_t min = 0.8;
1172     Float_t max = 1.2;
1173
1174     posX->SetMinimum(min);
1175     posX->SetMaximum(max);
1176     posY->SetMinimum(min);
1177     posY->SetMaximum(max);
1178     posZ->SetMinimum(min);
1179     posZ->SetMaximum(max);
1180
1181     posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1182
1183     posX->GetYaxis()->SetTitleOffset(1.7);
1184     posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1185     posY->GetYaxis()->SetTitleOffset(1.7);
1186     posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1187     posZ->GetYaxis()->SetTitleOffset(1.7);
1188     posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
1189
1190     posZ->GetXaxis()->SetRangeUser(0, 1);
1191
1192     TString canvasName;
1193     canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
1194
1195     TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1196     canvas->Divide(3, 1);
1197
1198     canvas->cd(1);
1199     InitPad();
1200     posX->DrawCopy();
1201
1202     canvas->cd(2);
1203     InitPad();
1204     posY->DrawCopy();
1205
1206     canvas->cd(3);
1207     InitPad();
1208     posZ->DrawCopy();
1209
1210     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1211     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1212   }
1213 }
1214
1215 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
1216 {
1217   TFile::Open(fileName);
1218   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1219   dNdEtaCorrection->LoadHistograms();
1220
1221   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1222   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1223
1224   gene->GetZaxis()->SetRangeUser(0.3, 10);
1225   meas->GetZaxis()->SetRangeUser(0.3, 10);
1226   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
1227   gene->GetZaxis()->SetRange(0, 0);
1228   meas->GetZaxis()->SetRange(0, 0);
1229
1230   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1231   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1232   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
1233   gene->GetYaxis()->SetRange(0, 0);
1234   meas->GetYaxis()->SetRange(0, 0);
1235
1236   gene->GetXaxis()->SetRangeUser(-10, 10);
1237   meas->GetXaxis()->SetRangeUser(-10, 10);
1238   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
1239   gene->GetXaxis()->SetRange(0, 0);
1240   meas->GetXaxis()->SetRange(0, 0);
1241 }
1242
1243 void Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1244 {
1245   gSystem->Load("libPWG0base");
1246
1247   Track2Particle2DCreatePlots(fileName);
1248
1249   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1250   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1251   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
1252
1253   /* this reads them from the file
1254   TH2* corrYX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx"));
1255   TH2* corrZX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx"));
1256   TH2* corrZY = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy"));*/
1257
1258   Prepare2DPlot(corrYX);
1259   Prepare2DPlot(corrZX);
1260   Prepare2DPlot(corrZY);
1261
1262   const char* title = "";
1263   corrYX->SetTitle(title);
1264   corrZX->SetTitle(title);
1265   corrZY->SetTitle(title);
1266
1267   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1268   canvas->Divide(3, 1);
1269
1270   canvas->cd(1);
1271   InitPadCOLZ();
1272   corrYX->Draw("COLZ");
1273
1274   canvas->cd(2);
1275   InitPadCOLZ();
1276   corrZX->Draw("COLZ");
1277
1278   canvas->cd(3);
1279   InitPadCOLZ();
1280   corrZY->Draw("COLZ");
1281
1282   canvas->SaveAs(Form("Track2Particle2D_%d.gif", gMax));
1283   canvas->SaveAs(Form("Track2Particle2D_%d.eps", gMax));
1284 }
1285
1286 void CompareTrack2Particle2D()
1287 {
1288   gSystem->Load("libPWG0base");
1289
1290   Track2Particle2DCreatePlots("correction_maponly-positive.root");
1291
1292   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("pos_yx"));
1293   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("pos_zx"));
1294   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("pos_zy"));
1295
1296   Track2Particle2DCreatePlots("correction_maponly-negative.root");
1297
1298   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("neg_yx"));
1299   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("neg_zx"));
1300   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("neg_zy"));
1301
1302   posYX->Divide(negYX);
1303   posZX->Divide(negZX);
1304   posZY->Divide(negZY);
1305
1306   Prepare2DPlot(posYX);
1307   Prepare2DPlot(posZX);
1308   Prepare2DPlot(posZY);
1309
1310   Float_t min = 0.8;
1311   Float_t max = 1.2;
1312
1313   posYX->SetMinimum(min);
1314   posYX->SetMaximum(max);
1315   posZX->SetMinimum(min);
1316   posZX->SetMaximum(max);
1317   posZY->SetMinimum(min);
1318   posZY->SetMaximum(max);
1319
1320   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1321   canvas->Divide(3, 1);
1322
1323   canvas->cd(1);
1324   InitPadCOLZ();
1325   posYX->Draw("COLZ");
1326
1327   canvas->cd(2);
1328   InitPadCOLZ();
1329   posZX->Draw("COLZ");
1330
1331   canvas->cd(3);
1332   InitPadCOLZ();
1333   posZY->Draw("COLZ");
1334
1335   canvas->SaveAs("CompareTrack2Particle2D.gif");
1336   canvas->SaveAs("CompareTrack2Particle2D.eps");
1337 }
1338
1339 void Track2Particle3D()
1340 {
1341   // get left margin proper
1342
1343   TFile* file = TFile::Open("correction_map.root");
1344
1345   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1346
1347   corr->SetTitle("Correction Factor");
1348   SetRanges(corr->GetZaxis());
1349
1350   Prepare3DPlot(corr);
1351
1352   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1353   canvas->SetTheta(29.428);
1354   canvas->SetPhi(16.5726);
1355
1356   corr->Draw();
1357
1358   canvas->SaveAs("Track2Particle3D.gif");
1359   canvas->SaveAs("Track2Particle3D.eps");
1360 }
1361
1362 void Track2Particle3DAll()
1363 {
1364   TFile* file = TFile::Open("correction_map.root");
1365
1366   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1367   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1368   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1369
1370   gene->SetTitle("Generated Particles");
1371   meas->SetTitle("Measured Tracks");
1372   corr->SetTitle("Correction Factor");
1373
1374   Prepare3DPlot(gene);
1375   Prepare3DPlot(meas);
1376   Prepare3DPlot(corr);
1377
1378   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1379   canvas->Divide(3, 1);
1380
1381   canvas->cd(1);
1382   InitPad();
1383   gene->Draw();
1384
1385   canvas->cd(2);
1386   meas->Draw();
1387
1388   canvas->cd(3);
1389   corr->Draw();
1390
1391   canvas->SaveAs("Track2Particle3DAll.gif");
1392   canvas->SaveAs("Track2Particle3DAll.eps");
1393 }
1394
1395 void MultiplicityMC(Int_t xRangeMax = 50)
1396 {
1397   TFile* file = TFile::Open("multiplicityMC.root");
1398
1399   if (!file)
1400   {
1401     printf("multiplicityMC.root could not be opened.\n");
1402     return;
1403   }
1404
1405   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1406   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1407   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1408
1409   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1410   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1411   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1412   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1413   {
1414     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1415     proj->Fit("gaus", "0");
1416     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1417     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1418
1419     continue;
1420
1421     // draw for debugging
1422     new TCanvas;
1423     proj->DrawCopy();
1424     proj->GetFunction("gaus")->DrawCopy("SAME");
1425   }
1426
1427   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1428
1429   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1430   {
1431     Float_t mean = correction->GetBinContent(i);
1432     Float_t width = correctionWidth->GetBinContent(i);
1433
1434     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1435     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1436     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1437
1438     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1439     {
1440       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1441     }
1442   }
1443
1444   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1445   fMultiplicityESDCorrectedRebinned->Rebin(10);
1446   fMultiplicityESDCorrectedRebinned->Scale(0.1);
1447
1448   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1449   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1450   ratio->Divide(fMultiplicityMC);
1451
1452   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1453   ratio2->Divide(fMultiplicityMC);
1454
1455   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1456   canvas->Divide(3, 2);
1457
1458   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1459   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1460   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1461   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1462   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1463   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1464   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1465
1466   canvas->cd(1); //InitPad();
1467   fMultiplicityESD->Draw();
1468   fMultiplicityMC->SetLineColor(2);
1469   fMultiplicityMC->Draw("SAME");
1470
1471   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1472   legend->AddEntry(fMultiplicityESD, "ESD");
1473   legend->AddEntry(fMultiplicityMC, "MC");
1474   legend->Draw();
1475
1476   canvas->cd(2);
1477   fCorrelation->Draw("COLZ");
1478
1479   canvas->cd(3);
1480   correction->Draw();
1481   //correction->Fit("pol1");
1482   correctionWidth->SetLineColor(2);
1483   correctionWidth->Draw("SAME");
1484
1485   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1486   legend->AddEntry(correction, "#bar{x}");
1487   legend->AddEntry(correctionWidth, "#sigma");
1488   legend->Draw();
1489
1490   canvas->cd(4);
1491   ratio->Draw();
1492
1493   ratio2->SetLineColor(2);
1494   ratio2->Draw("SAME");
1495
1496   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1497   legend->AddEntry(ratio, "uncorrected");
1498   legend->AddEntry(ratio2, "corrected");
1499   legend->Draw();
1500
1501   canvas->cd(5);
1502   fMultiplicityESDCorrected->SetLineColor(kBlue);
1503   fMultiplicityESDCorrected->Draw();
1504   fMultiplicityMC->Draw("SAME");
1505   fMultiplicityESD->Draw("SAME");
1506
1507   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1508   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1509   legend->AddEntry(fMultiplicityMC, "MC");
1510   legend->AddEntry(fMultiplicityESD, "ESD");
1511   legend->Draw();
1512
1513   canvas->cd(6);
1514   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
1515   fMultiplicityESDCorrectedRebinned->Draw();
1516   fMultiplicityMC->Draw("SAME");
1517
1518   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1519   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1520   legend->AddEntry(fMultiplicityMC, "MC");
1521   legend->Draw();
1522
1523   canvas->SaveAs("MultiplicityMC.gif");
1524 }
1525
1526 void MultiplicityESD()
1527 {
1528   TFile* file = TFile::Open("multiplicityESD.root");
1529
1530   if (!file)
1531   {
1532     printf("multiplicityESD.root could not be opened.\n");
1533     return;
1534   }
1535
1536   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1537
1538   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1539
1540   fMultiplicityESD->Draw();
1541 }
1542
1543 void drawPlots(Int_t max)
1544 {
1545   gMax = max;
1546
1547   ptCutoff();
1548   TriggerBias();
1549   VtxRecon();
1550   Track2Particle2D();
1551   Track2Particle3D();
1552   ptSpectrum();
1553   dNdEta();
1554 }
1555
1556 void drawPlots()
1557 {
1558   drawPlots(5);
1559   drawPlots(2);
1560 }
1561
1562 void CompareCorrection2Measured(const char* dataInput = "analysis_esd_raw.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1563 {
1564   loadlibs();
1565
1566   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1567   TFile::Open(correctionMapFile);
1568   dNdEtaCorrection->LoadHistograms();
1569
1570   TFile* file = TFile::Open(dataInput);
1571
1572   if (!file)
1573   {
1574     cout << "Error. File not found" << endl;
1575     return;
1576   }
1577
1578   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1579   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1580
1581   gROOT->cd();
1582
1583   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
1584   hist1->SetTitle("mc");
1585   Printf("mc contains %f entries", hist1->Integral());
1586   Printf("mc contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), 1, hist1->GetNbinsY(), hist1->GetZaxis()->FindBin(0.301), hist1->GetNbinsZ()));
1587
1588   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
1589   hist2->SetTitle("esd");
1590   Printf("esd contains %f entries", hist2->Integral());
1591   Printf("esd contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), 1, hist2->GetNbinsY(), hist2->GetZaxis()->FindBin(0.301), hist2->GetNbinsZ()));
1592
1593   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1594   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
1595
1596   hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1597   hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1598   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
1599
1600   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
1601   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
1602   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
1603   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
1604   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
1605 }
1606
1607 void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
1608 {
1609   loadlibs();
1610
1611   TFile* file = TFile::Open(dataInput);
1612
1613   if (!file)
1614   {
1615     cout << "Error. File not found" << endl;
1616     return;
1617   }
1618
1619   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1620   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1621
1622   TFile* file = TFile::Open(dataInput2);
1623
1624   if (!file)
1625   {
1626     cout << "Error. File not found" << endl;
1627     return;
1628   }
1629
1630   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
1631   fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
1632
1633   gROOT->cd();
1634
1635   TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
1636   hist1->SetTitle("esd1");
1637   Printf("esd1 contains %f entries", hist1->GetEntries());
1638   Printf("esd1 contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), 1, hist1->GetNbinsY(), hist1->GetZaxis()->FindBin(0.301), hist1->GetNbinsZ()));
1639
1640   TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
1641   hist2->SetTitle("esd2");
1642   Printf("esd2 contains %f entries", hist2->GetEntries());
1643   Printf("esd2 contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), 1, hist2->GetNbinsY(), hist2->GetZaxis()->FindBin(0.301), hist2->GetNbinsZ()));
1644
1645   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1646
1647   new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
1648   new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
1649   new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
1650
1651   TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
1652   TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
1653
1654   Printf("event1 contains %f entries", event1->GetEntries());
1655   Printf("event2 contains %f entries", event2->GetEntries());
1656   Printf("event1 integral is %f", event1->Integral());
1657   Printf("event2 integral is %f", event2->Integral());
1658   Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
1659   Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
1660
1661   projx1 = event1->ProjectionX();
1662   projx2 = event2->ProjectionX();
1663
1664   new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
1665
1666   projx1->Divide(projx2);
1667   new TCanvas; projx1->Draw();
1668
1669   event1->Divide(event2);
1670   new TCanvas; event1->Draw("COLZ");
1671
1672 }
1673
1674 void DrawTrackletOrigin()
1675 {
1676   TFile::Open("correction_map.root");
1677
1678   Int_t colors[]  = {1,2,3,4,6,7,8,102};
1679
1680   Int_t maxHists = 8;
1681   TH1* hist[8];
1682
1683   const char* titles[] = { "PP", "SS", "PP'", "PS", "PS*", "SP", "SS'", "" };
1684
1685   TLegend* legend = new TLegend(0.75, 0.6, 0.95, 0.95);
1686
1687   Int_t total = 0;
1688   for (Int_t i=0; i<maxHists; i++)
1689   {
1690     hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
1691     //hist[i]->Rebin(20);
1692     hist[i]->SetStats(kFALSE);
1693     hist[i]->SetLineColor(colors[i]);
1694     hist[i]->GetXaxis()->SetRangeUser(-0.2, 0.2);
1695     hist[i]->Draw(((i == 0) ? "" : "SAME"));
1696
1697     total += hist[i]->GetEntries();
1698
1699     if (i != 7)
1700       legend->AddEntry(hist[i], titles[i]);
1701   }
1702
1703   legend->Draw();
1704   gPad->SetLogy();
1705
1706   Printf("Total: %d", total);
1707   for (Int_t i=0; i<maxHists; i++)
1708     Printf("Histogram %d (%s) containts %.2f %% of the entries", i, titles[i], 100.0 * hist[i]->GetEntries() / total);
1709
1710   printf("|  Delta phi  |  Acc. %%  |  ");
1711   for (Int_t i=0; i<maxHists; i++)
1712     printf("%3s %%   |  ", titles[i]);
1713   Printf("");
1714
1715   for (Float_t f = 0.01; f < 0.09; f += 0.01)
1716   {
1717     Int_t integralBegin = hist[0]->GetXaxis()->FindBin(-f);
1718     Int_t integralEnd = hist[0]->GetXaxis()->FindBin(f);
1719
1720     Int_t total2 = 0;
1721     for (Int_t i=0; i<maxHists; i++)
1722       total2 += (Int_t) hist[i]->Integral(integralBegin, integralEnd);
1723
1724     printf("|    %.2f     |  %6.2f  |  ", f, 100.0 * total2 / total);
1725
1726     for (Int_t i=0; i<maxHists; i++)
1727       printf("%6.2f  |  ", (hist[i]->GetEntries() > 0) ? (100.0 * hist[i]->Integral(integralBegin, integralEnd) / hist[i]->GetEntries()) : -1.0);
1728     Printf("");
1729   }
1730 }
1731
1732 void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
1733 {
1734   loadlibs();
1735
1736   TFile::Open(fileName);
1737
1738   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
1739   if (!dNdEtaCorrection->LoadHistograms())
1740     return;
1741
1742   hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
1743
1744   proj = hist->Project3D("yx");
1745   proj->Draw("COLZ");
1746
1747   const Float_t limit = 10;
1748
1749   TString array = "{";
1750   TString arrayEnd = "}";
1751
1752   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
1753   {
1754     Int_t begin = -1;
1755     Int_t end = -1;
1756     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
1757     {
1758       if (begin == -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
1759         begin = x;
1760       if (begin != -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
1761         end = x;
1762     }
1763     Printf("Limits for y = %d are %d to %d", y, begin, end);
1764
1765     if (y > 1)
1766       array += ", ";
1767     array += Form("%d", begin);
1768
1769     if (y > 1)
1770       arrayEnd.Prepend(", ");
1771     arrayEnd.Prepend(Form("%d", (end == -1) ? -1 : proj->GetNbinsX() + 1 - end));
1772   }
1773   array += "}";
1774   arrayEnd.Prepend("{");
1775
1776   Printf("Begin array:");
1777   Printf("%s", array.Data());
1778
1779   Printf("End array (mirrored) (should be the same):");
1780   Printf("%s", arrayEnd.Data());
1781 }
1782
1783 void AverageMultiplicity(const char* fileName = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1784 {
1785   loadlibs();
1786
1787   TFile::Open(fileName);
1788
1789   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1790   dNdEtaCorrection->LoadHistograms();
1791   TH2* events = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram();
1792   TH3* tracks = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram();
1793
1794   Float_t nEvents = events->Integral(events->GetXaxis()->FindBin(-1), events->GetXaxis()->FindBin(1), 0, events->GetNbinsY()+1);
1795   Float_t nTracks = tracks->Integral(tracks->GetXaxis()->FindBin(-1), tracks->GetXaxis()->FindBin(1), tracks->GetYaxis()->FindBin(-0.39), tracks->GetYaxis()->FindBin(0.59), 0, tracks->GetNbinsZ()+1);
1796
1797   Printf("%f %f --> %f", nEvents, nTracks, nTracks / nEvents);
1798 }
1799