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("libVMC");
  26   gSystem->Load("libTree");
  27   gSystem->Load("libSTEERBase");
  28   gSystem->Load("libESD");
  29   gSystem->Load("libAOD");
  30   gSystem->Load("libANALYSIS");
  31   gSystem->Load("libANALYSISalice");
  32   gSystem->Load("libPWG0base");
  33 }
  34
  35 void SetRanges(TAxis* axis)
  36 {
  37   if (strcmp(axis->GetTitle(), "#eta") == 0)
  38     axis->SetRangeUser(-1.4999, 1.4999);
  39   if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0 || strcmp(axis->GetTitle(), "p_{T} (GeV/c)") == 0)
  40   {
  41     axis->SetRangeUser(0, 4.9999);
  42     axis->SetTitle("p_{T} (GeV/c)");
  43   }
  44   if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0 || strcmp(axis->GetTitle(), "vtx z (cm)") == 0)
  45   {
  46     axis->SetRangeUser(-15, 14.9999);
  47     axis->SetTitle("vtx-z (cm)");
  48   }
  49   if (strcmp(axis->GetTitle(), "Ntracks") == 0)
  50     axis->SetRangeUser(0, 99.9999);
  51 }
  52
  53 void SetRanges(TH1* hist)
  54 {
  55   SetRanges(hist->GetXaxis());
  56   SetRanges(hist->GetYaxis());
  57   SetRanges(hist->GetZaxis());
  58 }
  59
  60
  61 void Prepare3DPlot(TH3* hist)
  62 {
  63   hist->GetXaxis()->SetTitleOffset(1.5);
  64   hist->GetYaxis()->SetTitleOffset(1.5);
  65   hist->GetZaxis()->SetTitleOffset(1.5);
  66
  67   hist->SetStats(kFALSE);
  68 }
  69
  70 void Prepare2DPlot(TH2* hist)
  71 {
  72   hist->SetStats(kFALSE);
  73   hist->GetYaxis()->SetTitleOffset(1.4);
  74
  75   hist->SetMinimum(0);
  76   hist->SetMaximum(gMax);
  77
  78   SetRanges(hist);
  79 }
  80
  81 void Prepare1DPlot(TH1* hist)
  82 {
  83   hist->SetLineWidth(2);
  84   hist->SetStats(kFALSE);
  85
  86   hist->GetXaxis()->SetLabelOffset(0.02);
  87   hist->GetXaxis()->SetTitleOffset(1.3);
  88   hist->GetYaxis()->SetTitleOffset(1.3);
  89
  90   SetRanges(hist);
  91 }
  92
  93 void InitPad()
  94 {
  95   gPad->Range(0, 0, 1, 1);
  96   gPad->SetLeftMargin(0.15);
  97   //gPad->SetRightMargin(0.05);
  98   //gPad->SetTopMargin(0.13);
  99   gPad->SetBottomMargin(0.12);
 100
 101   gPad->SetGridx();
 102   gPad->SetGridy();
 103 }
 104
 105 void InitPadCOLZ()
 106 {
 107   gPad->Range(0, 0, 1, 1);
 108   gPad->SetRightMargin(0.15);
 109   gPad->SetLeftMargin(0.12);
 110   gPad->SetTopMargin(0.05);
 111
 112   gPad->SetGridx();
 113   gPad->SetGridy();
 114 }
 115
 116 // --- end of helpers --- begin functions ---
 117
 118 void DrawOverview(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
 119 {
 120   loadlibs();
 121   if (!TFile::Open(fileName))
 122     return;
 123
 124   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
 125   if (!dNdEtaCorrection->LoadHistograms())
 126     return;
 127
 128   dNdEtaCorrection->Finish();
 129
 130   dNdEtaCorrection->DrawOverview();
 131 }
 132
 133 void PrintInfo(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
 134 {
 135   loadlibs();
 136   if (!TFile::Open(fileName))
 137     return;
 138
 139   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
 140   if (!dNdEtaCorrection->LoadHistograms())
 141     return;
 142
 143   dNdEtaCorrection->Finish();
 144
 145   for (Int_t i=AlidNdEtaCorrection::kTrack2Particle; i<=AlidNdEtaCorrection::kND; i++)
 146   {
 147     Printf("Correction %d", i);
 148     dNdEtaCorrection->GetCorrection(i)->PrintInfo(0.3);
 149   }
 150 }
 151
 152 void PrintAllInfos()
 153 {
 154   PrintInfo();
 155
 156   Printf("RAW ESD");
 157   TFile::Open("analysis_esd_raw.root");
 158   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 159   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
 160   fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
 161
 162   const Int_t num = 3;
 163   const char* results[] = { "dndeta", "dndetaTr", "dndetaTrVtx" };
 164
 165   TFile::Open("analysis_esd.root");
 166   for (Int_t i=0; i<num; i++)
 167   {
 168     Printf("CORRECTED %s", results[i]);
 169     dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 170     fdNdEtaAnalysis->LoadHistograms(results[i]);
 171     fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
 172   }
 173 }
 174
 175 void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
 176 {
 177   gSystem->Load("libPWG0base");
 178
 179   TFile::Open(esdFile);
 180   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 181   fdNdEtaAnalysisESD->LoadHistograms();
 182
 183   TFile::Open(mcFile);
 184   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
 185   fdNdEtaAnalysisMC->LoadHistograms();
 186   //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 187
 188   for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
 189     fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
 190
 191   fdNdEtaAnalysisESD->DrawHistograms();
 192 }
 193
 194 void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 195 {
 196   gSystem->Load("libPWG0base");
 197
 198   const char* ESDfolder = 0;
 199
 200   if (plot == 0) // all
 201     ESDfolder = "dndeta";
 202   else if (plot == 1) // mb
 203     ESDfolder = "dndeta_mb";
 204   else if (plot == 2) // mb vtx
 205     ESDfolder = "dndeta_mbvtx";
 206
 207   TFile::Open("analysis_esd.root");
 208   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 209   fdNdEtaAnalysisESD->LoadHistograms();
 210
 211   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 212   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 213
 214   TH2F* hist = 0;
 215
 216   if (plot == 0) // all
 217     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
 218   else if (plot == 1) // mb
 219     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
 220   else if (plot == 2) // mb vtx
 221     hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
 222
 223   TH1* proj = hist->ProjectionX();
 224
 225   TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
 226   for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
 227   {
 228     Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
 229     if (value != 0)
 230     {
 231       printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
 232       vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
 233     }
 234   }
 235
 236   new TCanvas;
 237   vertex->DrawCopy();
 238 }
 239
 240 void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 241 {
 242   gSystem->Load("libPWG0base");
 243
 244   TFile::Open("analysis_esd.root");
 245   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
 246   fdNdEtaAnalysisESD->LoadHistograms();
 247
 248   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 249   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 250
 251   //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
 252   //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 253
 254   TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
 255   TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 256
 257   TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
 258
 259   new TCanvas;
 260   histESD->Draw();
 261
 262   new TCanvas;
 263   histCorr->Draw();
 264
 265   for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
 266     for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
 267       for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
 268       {
 269         Float_t value1 = histESD->GetBinContent(x, y, z);
 270         Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
 271
 272         if (value2 > 0 && value1 > 0)
 273         {
 274           printf("%f %f %f\n", value1, value2, value1 / value2);
 275           diff->Fill(value1 / value2);
 276         }
 277       }
 278
 279   new TCanvas;
 280   diff->Draw();
 281 }
 282
 283 Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
 284 {
 285   Double_t avgMC = 0;
 286   Double_t avgESD = 0;
 287   for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
 288   {
 289     avgMC += histMC->GetBinContent(bin);
 290     avgESD += histESD->GetBinContent(bin);
 291   }
 292   Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
 293
 294   avgMC /= nBins;
 295   avgESD /= nBins;
 296
 297   // deviation when integrate in |eta| < 0.8 between mc and esd
 298   Double_t diffFullRange = (avgMC - avgESD) / avgMC;
 299
 300   Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
 301
 302   return diffFullRange;
 303 }
 304
 305 void dNdEtaNoResolution()
 306 {
 307   loadlibs();
 308
 309   TFile::Open("correction_map.root");
 310
 311   const char* correctionMapFolder = "dndeta_correction";
 312   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 313   dNdEtaCorrection->LoadHistograms();
 314   dNdEtaCorrection->GetTrack2ParticleCorrection()->PrintInfo(0.3);
 315
 316   TFile::Open("analysis_mc.root");
 317   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
 318   fdNdEtaAnalysis->LoadHistograms();
 319   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, "ESD (no resolution effect) -> MB with vertex");
 320   fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0)->SetMarkerStyle(21);
 321
 322   TFile::Open("analysis_mc.root");
 323   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
 324   fdNdEtaAnalysisMC->LoadHistograms();
 325   fdNdEtaAnalysisMC->Finish(0, 0, AlidNdEtaCorrection::kNone, "MC: MB with vertex");
 326
 327   DrawdNdEtaRatio(fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0), fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(0), "MB with vertex (no resolution effect)", 3);
 328 }
 329
 330 TH1* GetMCHist(const char* folder, Float_t ptCut, const char* tag)
 331 {
 332   loadlibs();
 333
 334   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(folder, folder);
 335   fdNdEtaAnalysis->LoadHistograms();
 336   fdNdEtaAnalysis->Finish(0, ptCut, AlidNdEtaCorrection::kNone, tag);
 337   return fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 338 }
 339
 340 void dNdEtaFinal(Bool_t spd = kTRUE)
 341 {
 342   TFile* file = TFile::Open("analysis_esd.root");
 343   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
 344   TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
 345   Prepare1DPlot(histESD);
 346   Prepare1DPlot(histESDnsd);
 347
 348   TCanvas* canvas = new TCanvas("dNdEtaFinal", "dNdEtaFinal", 600, 600);
 349   gPad->SetTopMargin(0.05);
 350   gPad->SetRightMargin(0.05);
 351   gPad->SetLeftMargin(0.12);
 352   gPad->SetBottomMargin(0.12);
 353   gPad->SetGridx();
 354   gPad->SetGridy();
 355
 356   Float_t etaMax = 1.9;
 357   Float_t histMax = 1.39;
 358   Float_t systErrorValue = 0.023;
 359   Float_t systErrorNSDValue = 0.081;
 360   if (!spd)
 361   {
 362     //etaMax = 1.5;
 363     histMax = 0.99;
 364     systErrorValue = 0.043;
 365     systErrorNSDValue = 0.088;
 366   }
 367
 368   dummy = new TH2F("dummy", ";#eta;dN_{ch}/d#eta", 100, -etaMax, etaMax, 100, 3, 8);
 369   dummy->SetStats(0);
 370   dummy->GetYaxis()->SetTitleOffset(1.3);
 371
 372   histESD->SetMarkerStyle(20);
 373   histESDnsd->SetMarkerStyle(21);
 374   histESDnsd->SetMarkerColor(4);
 375   histESDnsd->SetLineColor(4);
 376   histESD->SetMarkerSize(1.5);
 377   histESDnsd->SetMarkerSize(1.5);
 378
 379   histESD->GetXaxis()->SetRangeUser(-histMax, histMax);
 380   histESDnsd->GetXaxis()->SetRangeUser(-histMax, histMax);
 381
 382   legend = new TLegend(0.3, 0.2, 0.78, 0.4);
 383   legend->SetFillColor(0);
 384   legend->SetTextSize(0.04);
 385   legend->AddEntry(histESD, "Inelastic events", "P");
 386   legend->AddEntry(histESDnsd, "NSD events", "P");
 387
 388   dummy->Draw();
 389
 390   // syst errors.
 391   TH1* systError = (TH1*) histESD->Clone("systError");
 392   for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
 393     systError->SetBinError(i, systError->GetBinContent(i) * systErrorValue);
 394   // change error drawing style
 395   systError->SetFillColor(15);
 396   systError->DrawCopy("SAME E2 ][");
 397
 398   // syst error NSD
 399   for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
 400   {
 401     systError->SetBinContent(i, histESDnsd->GetBinContent(i));
 402     systError->SetBinError(i, systError->GetBinContent(i) * systErrorNSDValue);
 403   }
 404   systError->DrawCopy("SAME E2 ][");
 405
 406   histESD->Draw("SAME");
 407   histESDnsd->Draw("SAME");
 408   legend->Draw();
 409
 410   canvas->SaveAs(Form("%s_dndeta_final.eps", (spd) ? "spd" : "tpc"));
 411 }
 412
 413 void dNdEtaPythiaPhojet()
 414 {
 415   // evtl. deactivate acceptance maps in dNdEtaAnalysis.cxx
 416
 417   loadlibs();
 418
 419   TH1* hist[4];
 420
 421   TFile::Open("LHC08c11_10TeV_0.5T/mb1/spd/analysis_mc.root");
 422   hist[0] =         (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
 423   hist[1] =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
 424
 425   TFile::Open("LHC08c15_10TeV_0.5T_Phojet/mb1/spd/analysis_mc.root");
 426   hist[2] =         (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMCPhojet");
 427   hist[3] =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsdPhojet");
 428
 429   file = TFile::Open("pythia_phojet_dndeta.root", "RECREATE");
 430   for (Int_t i=0; i<4; i++)
 431     hist[i]->Write();
 432   file->Close();
 433 }
 434
 435 void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 436 {
 437   loadlibs();
 438
 439   TFile* file = TFile::Open("analysis_esd.root");
 440
 441   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 442   fdNdEtaAnalysis->LoadHistograms("dndeta");
 443
 444   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
 445   TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
 446   TH1* histESDnsdNoPt = (TH1*) file->Get("dndetaNSD/dNdEta");
 447   TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
 448   TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
 449   TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
 450   TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
 451   TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
 452   TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
 453
 454   Prepare1DPlot(histESD);
 455   Prepare1DPlot(histESDnsd);
 456   Prepare1DPlot(histESDMB);
 457   Prepare1DPlot(histESDMBVtx);
 458
 459   Prepare1DPlot(histESDNoPt);
 460   Prepare1DPlot(histESDMBNoPt);
 461   Prepare1DPlot(histESDMBVtxNoPt);
 462   Prepare1DPlot(histESDMBTracksNoPt);
 463
 464   histESD->SetLineWidth(0);
 465   histESDnsd->SetLineWidth(0);
 466   histESDMB->SetLineWidth(0);
 467   histESDMBVtx->SetLineWidth(0);
 468
 469   histESDNoPt->SetLineWidth(0);
 470   histESDMBNoPt->SetLineWidth(0);
 471   histESDMBVtxNoPt->SetLineWidth(0);
 472
 473   histESD->SetMarkerColor(1);
 474   histESDnsd->SetMarkerColor(6);
 475   histESDMB->SetMarkerColor(2);
 476   histESDMBVtx->SetMarkerColor(4);
 477
 478   histESD->SetLineColor(1);
 479   histESDnsd->SetLineColor(6);
 480   histESDMB->SetLineColor(2);
 481   histESDMBVtx->SetLineColor(4);
 482
 483   histESDNoPt->SetMarkerColor(1);
 484   histESDMBNoPt->SetMarkerColor(2);
 485   histESDMBVtxNoPt->SetMarkerColor(4);
 486   histESDMBTracksNoPt->SetMarkerColor(3);
 487
 488   histESD->SetMarkerStyle(20);
 489   histESDnsd->SetMarkerStyle(29);
 490   histESDMB->SetMarkerStyle(21);
 491   histESDMBVtx->SetMarkerStyle(22);
 492
 493   histESDNoPt->SetMarkerStyle(20);
 494   histESDMBNoPt->SetMarkerStyle(21);
 495   histESDMBVtxNoPt->SetMarkerStyle(22);
 496   histESDMBTracksNoPt->SetMarkerStyle(23);
 497
 498   Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 0.89 : 1.79;
 499   Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 1.2 : 2.3;
 500   //Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 0.89 : 1.39;
 501   //Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 1.2 : 1.9;
 502
 503   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 504   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 505   histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 506   histESDnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 507
 508   histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 509   histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 510   histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 511   histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 512
 513   Float_t max = TMath::Max(histESDMBVtx->GetMaximum(), histESDMB->GetMaximum());
 514   max = TMath::Max(max, histESD->GetMaximum());
 515
 516   TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.4);
 517   legend->SetFillColor(0);
 518   legend->AddEntry(histESDMBVtx, "Triggered, vertex");
 519   legend->AddEntry(histESDMB, "Triggered");
 520   legend->AddEntry(histESD, "All events");
 521
 522   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 2.1, max * 1.1);
 523   Prepare1DPlot(dummy);
 524   dummy->SetStats(kFALSE);
 525   dummy->SetXTitle("#eta");
 526   dummy->SetYTitle("dN_{ch}/d#eta");
 527   dummy->GetYaxis()->SetTitleOffset(1);
 528
 529   TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
 530
 531   dummy->DrawCopy();
 532   histESDMBVtx->Draw("SAME");
 533   histESDMB->Draw("SAME");
 534   histESD->Draw("SAME");
 535   legend->Draw();
 536
 537   if (save)
 538   {
 539     canvas->SaveAs("dNdEta1.gif");
 540     canvas->SaveAs("dNdEta1.eps");
 541   }
 542
 543   if (onlyESD)
 544     return;
 545
 546   loadlibs();
 547
 548   TFile* file2 = TFile::Open("analysis_mc.root");
 549
 550   TH1* histMCTrVtx =       (TH1*) GetMCHist("dndetaTrVtx", -1, "MC: MB with trigger")->Clone("histMCTrVtx");
 551   TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
 552
 553   TH1* histMC =            (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
 554   TH1* histMCTr =          (TH1*) GetMCHist("dndetaTr", -1, "MC: minimum bias")->Clone("histMCTr");
 555   TH1* histMCnsd =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
 556
 557   TH1* histMCPtCut =       (TH1*) GetMCHist("dndeta", 0.21, "MC: full inelastic, pt cut")->Clone("histMCPtCut");
 558   TH1* histMCTrPtCut =     (TH1*) GetMCHist("dndetaTr", 0.21, "MC: minimum bias, pt cut")->Clone("histMCTrPtCut");
 559   TH1* histMCTrVtxPtCut =  (TH1*) GetMCHist("dndetaTrVtx", 0.21, "MC: MB with trigger, pt cut")->Clone("histMCTrVtxPtCut");
 560   TH1* histMCnsdNoPt =     (TH1*) GetMCHist("dndetaNSD", 0.21, "MC: NSD, put cut")->Clone("histMCnsdNoPt");
 561   TH1* histMCTracksPtCut = (TH1*) GetMCHist("dndetaTracks", 0.21, "MC: Tracks w/o resolution effect, pt cut")->Clone("histMCTracksPtCut");
 562
 563   Prepare1DPlot(histMC);
 564   Prepare1DPlot(histMCnsd);
 565   Prepare1DPlot(histMCTr);
 566   Prepare1DPlot(histMCTrVtx);
 567
 568   Prepare1DPlot(histMCPtCut);
 569   Prepare1DPlot(histMCTrPtCut);
 570   Prepare1DPlot(histMCTrVtxPtCut);
 571   Prepare1DPlot(histMCTracksPtCut);
 572
 573   histMC->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 574   histMCnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 575   histMCTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 576   histMCTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 577
 578   histMCPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 579   histMCTrPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 580   histMCTrVtxPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 581   histMCTracksPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 582
 583   histMC->SetLineColor(1);
 584   histMCnsd->SetLineColor(6);
 585   histMCTr->SetLineColor(2);
 586   histMCTrVtx->SetLineColor(4);
 587
 588   histMCPtCut->SetLineColor(1);
 589   histMCTrPtCut->SetLineColor(2);
 590   histMCTrVtxPtCut->SetLineColor(4);
 591   if (histMCTracksPtCut)
 592     histMCTracksPtCut->SetLineColor(3);
 593
 594   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
 595
 596   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
 597   dummy2->GetYaxis()->SetRangeUser(0, max * 1.1);
 598
 599   dummy2->DrawCopy();
 600   histMC->Draw("SAME");
 601   histMCnsd->Draw("SAME");
 602   histMCTr->Draw("SAME");
 603   histMCTrVtx->Draw("SAME");
 604   histESD->Draw("SAME");
 605   histESDnsd->Draw("SAME");
 606   histESDMB->Draw("SAME");
 607   histESDMBVtx->Draw("SAME");
 608   histESDNoPt->Draw("SAME");
 609   histESDMBNoPt->Draw("SAME");
 610   histESDMBVtxNoPt->Draw("SAME");
 611   histESDMBTracksNoPt->Draw("SAME");
 612   histMCPtCut->Draw("SAME");
 613   histMCTrPtCut->Draw("SAME");
 614   histMCTrVtxPtCut->Draw("SAME");
 615   if (histMCTracksPtCut)
 616     histMCTracksPtCut->Draw("SAME");
 617
 618   if (save)
 619   {
 620     canvas2->SaveAs("dNdEta2.gif");
 621     canvas2->SaveAs("dNdEta2.eps");
 622   }
 623
 624   TH1* ratio = (TH1*) DrawdNdEtaRatio(histESD, histMC, "full_inelastic", etaPlotLimit)->Clone();
 625   TH1* ratioTr = (TH1*) DrawdNdEtaRatio(histESDMB, histMCTr, "triggered", etaPlotLimit)->Clone();
 626   TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
 627   TH1* ratioTrVtxNoPt = (TH1*) DrawdNdEtaRatio(histESDMBVtxNoPt, histMCTrVtxPtCut, "triggered_vertex_nopt", etaPlotLimit)->Clone();
 628   TH1* ratioNSD = (TH1*) DrawdNdEtaRatio(histESDnsd, histMCnsd, "NSD", etaPlotLimit)->Clone();
 629
 630   // draw ratios of single steps
 631   c7 = new TCanvas("all_ratios", "all_ratios", 600, 600);
 632   c7->SetRightMargin(0.05);
 633   c7->SetTopMargin(0.05);
 634   c7->SetGridx();
 635   c7->SetGridy();
 636
 637   ratioTrVtxNoPt->SetMarkerStyle(20);
 638   ratioTrVtx->SetMarkerStyle(21);
 639   ratioTr->SetMarkerStyle(23);
 640   ratio->SetMarkerStyle(22);
 641   ratioNSD->SetMarkerStyle(26);
 642
 643   ratioTrVtxNoPt->SetMarkerSize(2);
 644   ratioTrVtx->SetMarkerSize(2);
 645   ratioTr->SetMarkerSize(2);
 646   ratio->SetMarkerSize(2);
 647   ratioNSD->SetMarkerSize(2);
 648
 649   ratioTrVtxNoPt->SetMarkerColor(1);
 650   ratioTrVtx->SetMarkerColor(2);
 651   ratioTr->SetMarkerColor(4);
 652   ratio->SetMarkerColor(2);
 653   ratioNSD->SetMarkerColor(1);
 654
 655   ratioTrVtxNoPt->SetLineColor(1);
 656   ratioTrVtx->SetLineColor(2);
 657   ratioTr->SetLineColor(4);
 658   ratio->SetLineColor(2);
 659   ratioNSD->SetLineColor(1);
 660
 661   legend7 = new TLegend(0.13, 0.7, 0.94, 0.9);
 662   legend7->SetFillColor(0);
 663   legend7->SetTextSize(0.035);
 664   legend7->SetNColumns(2);
 665
 666   flat = new TF1("flat", "-1", -5, 5);
 667   ratioTrVtxNoPt->Add(flat);
 668   ratioTrVtx->Add(flat);
 669   ratioTr->Add(flat);
 670   ratio->Add(flat);
 671   ratioNSD->Add(flat);
 672
 673   ratioTrVtxNoPt->Scale(100);
 674   ratioTrVtx->Scale(100);
 675   ratioTr->Scale(100);
 676   ratio->Scale(100);
 677   ratioNSD->Scale(100);
 678
 679   ratio->Add(ratioTr, -1);
 680   ratioNSD->Add(ratioTr, -1);
 681   ratioTr->Add(ratioTrVtx, -1);
 682   ratioTrVtx->Add(ratioTrVtxNoPt, -1);
 683
 684   legend7->AddEntry(ratioTrVtxNoPt, "Track-to-particle", "P");
 685   legend7->AddEntry(ratio, "Trigger-bias INEL", "P");
 686   legend7->AddEntry(ratioTr, "Vertex-reconstruction", "P");
 687   legend7->AddEntry(ratioNSD, "Trigger-bias NSD", "P");
 688   if (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC)
 689     legend7->AddEntry(ratioTrVtx, "p_{T} cut-off", "P");
 690
 691   TH1* dummy7 = new TH2F("dummy7", ";#eta;Deviation in %", 100, -etaPlotLimit, etaPlotLimit, 100, -5, 7);
 692   dummy7->SetStats(0);
 693   dummy7->Draw();
 694
 695   ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 696   ratioTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 697   ratioTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 698   ratioTrVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 699   ratioNSD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 700
 701   ratio->Draw("HIST EP SAME");
 702   ratioTr->Draw("HIST EP SAME");
 703   if (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC)
 704     ratioTrVtx->Draw("HIST EP SAME");
 705   ratioTrVtxNoPt->Draw("HIST EP SAME");
 706   ratioNSD->Draw("HIST EP SAME");
 707   legend7->Draw();
 708
 709   c7->SaveAs("ratios.eps");
 710
 711   new TCanvas;
 712   dummy2->DrawCopy();
 713   histMCnsd->Draw("SAME");
 714   histESDnsd->Draw("SAME");
 715
 716   ratio = (TH1*) histMC->Clone("ratio");
 717   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
 718
 719   ratio->Divide(histESD);
 720   ratioNoPt->Divide(histESDNoPt);
 721
 722   ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 723
 724   ratio->SetLineColor(1);
 725   ratioNoPt->SetLineColor(2);
 726
 727   Double_t average = 0;       // average deviation from 1 in ratio (depends on the number of bins if statistical)
 728   for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
 729     average += TMath::Abs(ratio->GetBinContent(bin) - 1);
 730   Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
 731   average /= nBins;
 732   Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
 733
 734   PrintIntegratedDeviation(histMC, histESD, "all events");
 735   PrintIntegratedDeviation(histMCnsd, histESDnsd, "all events (NSD)");
 736   PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
 737   PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
 738   PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
 739   PrintIntegratedDeviation(histMCnsdNoPt, histESDnsdNoPt, "all events (NSD) (no pt corr)");
 740   PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
 741   PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
 742
 743   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 600, 600);
 744   canvas3->Range(0, 0, 1, 1);
 745   //canvas3->Divide(1, 2, 0, 0);
 746
 747   //canvas3->cd(1);
 748   TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
 749   pad1->SetTopMargin(0.05);
 750   pad1->SetLeftMargin(0.13);
 751   pad1->Draw();
 752
 753   TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
 754   pad2->SetLeftMargin(0.13);
 755   pad2->Draw();
 756
 757   pad1->SetRightMargin(0.01);
 758   pad2->SetRightMargin(0.01);
 759
 760   // no border between them
 761   pad1->SetBottomMargin(0);
 762   pad2->SetTopMargin(0);
 763
 764   pad1->cd();
 765   pad1->SetGridx();
 766   pad1->SetGridy();
 767
 768   legend->AddEntry(histMC, "MC prediction");
 769
 770   dummy->GetXaxis()->SetLabelSize(0.08);
 771   dummy->GetYaxis()->SetLabelSize(0.08);
 772   dummy->GetXaxis()->SetTitleSize(0.08);
 773   dummy->GetYaxis()->SetTitleSize(0.08);
 774   dummy->GetYaxis()->SetTitleOffset(0.8);
 775   dummy->DrawCopy();
 776   histESDMBVtx->Draw("SAME");
 777   histESDMB->Draw("SAME");
 778   histESD->Draw("SAME");
 779   histMC->Draw("SAME");
 780
 781   legend->SetTextSize(0.08);
 782   legend->Draw();
 783
 784   pad2->cd();
 785   pad2->SetBottomMargin(0.15);
 786   //pad2->SetGridx();
 787   //pad2->SetGridy();
 788
 789   Float_t minR = 0.91; //TMath::Min(0.961, ratio->GetMinimum() * 0.95);
 790   Float_t maxR = 1.09; //TMath::Max(1.049, ratio->GetMaximum() * 1.05);
 791
 792   TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
 793   dummy3.SetStats(kFALSE);
 794   for (Int_t i=1; i<=100; ++i)
 795     dummy3.SetBinContent(i, 1);
 796   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
 797   dummy3.SetLineWidth(2);
 798   dummy3.GetXaxis()->SetLabelSize(0.08);
 799   dummy3.GetYaxis()->SetLabelSize(0.08);
 800   dummy3.GetXaxis()->SetTitleSize(0.08);
 801   dummy3.GetYaxis()->SetTitleSize(0.08);
 802   dummy3.GetYaxis()->SetTitleOffset(0.8);
 803   dummy3.DrawCopy();
 804
 805   ratio->Draw("SAME");
 806
 807   //pad2->Draw();
 808
 809   canvas3->Modified();
 810
 811   if (save)
 812   {
 813     canvas3->SaveAs("dNdEta.gif");
 814     canvas3->SaveAs("dNdEta.eps");
 815   }
 816
 817   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
 818
 819   ratio->Draw();
 820   ratioNoPt->Draw("SAME");
 821
 822   TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
 823   legend->SetFillColor(0);
 824   legend->AddEntry(ratio, "mc/esd");
 825   legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
 826   legend->Draw();
 827 }
 828
 829 TH1* DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
 830 {
 831   TCanvas* canvas3 = new TCanvas(name, name, 600, 600);
 832   canvas3->Range(0, 0, 1, 1);
 833
 834   TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
 835   pad1->Draw();
 836
 837   TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
 838   pad2->Draw();
 839
 840   pad1->SetRightMargin(0.01);
 841   pad2->SetRightMargin(0.01);
 842   pad1->SetTopMargin(0.05);
 843   pad1->SetLeftMargin(0.13);
 844   pad2->SetLeftMargin(0.13);
 845   pad2->SetBottomMargin(0.15);
 846
 847   // no border between them
 848   pad1->SetBottomMargin(0);
 849   pad2->SetTopMargin(0);
 850
 851   pad1->cd();
 852   pad1->SetGridx();
 853   pad1->SetGridy();
 854
 855   TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.3);
 856   legend->SetFillColor(0);
 857   legend->AddEntry(corr, "Corrected");
 858   legend->AddEntry(mc, "MC prediction");
 859   legend->SetTextSize(0.08);
 860
 861   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 3.1, corr->GetMaximum() * 1.1);
 862   Prepare1DPlot(dummy);
 863   dummy->SetStats(kFALSE);
 864   dummy->SetXTitle("#eta");
 865   dummy->SetYTitle("dN_{ch}/d#eta");
 866   dummy->GetYaxis()->SetTitleOffset(1);
 867
 868   dummy->GetXaxis()->SetLabelSize(0.08);
 869   dummy->GetYaxis()->SetLabelSize(0.08);
 870   dummy->GetXaxis()->SetTitleSize(0.08);
 871   dummy->GetYaxis()->SetTitleSize(0.08);
 872   dummy->GetYaxis()->SetTitleOffset(0.8);
 873   dummy->DrawCopy();
 874
 875   corr->Draw("SAME");
 876   mc->Draw("SAME");
 877
 878   legend->Draw();
 879
 880   pad2->cd();
 881   pad2->SetBottomMargin(0.15);
 882   //pad2->SetGridx();
 883   //pad2->SetGridy();
 884
 885   TH1* ratio = (TH1*) mc->Clone("ratio");
 886   ratio->Divide(corr);
 887
 888   Float_t minR = TMath::Min(0.91, ratio->GetMinimum() * 0.95);
 889   Float_t maxR = TMath::Max(1.09, ratio->GetMaximum() * 1.05);
 890
 891   TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
 892   dummy3.SetStats(kFALSE);
 893   for (Int_t i=1; i<=100; ++i)
 894         dummy3.SetBinContent(i, 1);
 895   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
 896   dummy3.SetLineWidth(2);
 897   dummy3.GetXaxis()->SetLabelSize(0.08);
 898   dummy3.GetYaxis()->SetLabelSize(0.08);
 899   dummy3.GetXaxis()->SetTitleSize(0.08);
 900   dummy3.GetYaxis()->SetTitleSize(0.08);
 901   dummy3.GetYaxis()->SetTitleOffset(0.8);
 902   dummy3.DrawCopy();
 903
 904   ratio->Draw("SAME");
 905
 906   canvas3->Modified();
 907
 908   return ratio;
 909 }
 910
 911 void ptSpectrum()
 912 {
 913   TFile* file = TFile::Open("analysis_esd.root");
 914   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
 915
 916   TFile* file2 = TFile::Open("analysis_mc.root");
 917   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
 918
 919   TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
 920   InitPad();
 921   gPad->SetLogy();
 922
 923   Prepare1DPlot(histMC);
 924   Prepare1DPlot(histESD);
 925
 926   histESD->SetTitle("");
 927   histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 928   histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
 929
 930   histMC->SetLineColor(kBlue);
 931   histESD->SetLineColor(kRed);
 932
 933   histESD->GetYaxis()->SetTitleOffset(1.5);
 934   histESD->GetXaxis()->SetRangeUser(0, 4.9999);
 935
 936   histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
 937
 938   histESD->Draw();
 939   histMC->Draw("SAME");
 940
 941   canvas->SaveAs("ptSpectrum.gif");
 942   canvas->SaveAs("ptSpectrum.eps");
 943 }
 944
 945 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
 946 {
 947   gSystem->Load("libPWG0base");
 948
 949   TFile::Open(fileName);
 950   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 951   dNdEtaCorrection->LoadHistograms();
 952
 953   TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
 954   TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
 955
 956   Prepare2DPlot(corrTrigger);
 957   corrTrigger->SetTitle("b) Trigger bias correction");
 958
 959   Prepare2DPlot(corrVtx);
 960   corrVtx->SetTitle("a) Vertex reconstruction correction");
 961
 962   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
 963   corrVtx->GetYaxis()->SetTitle("Multiplicity");
 964
 965   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
 966   canvas->Divide(2, 1);
 967
 968   canvas->cd(1);
 969   InitPadCOLZ();
 970   corrVtx->DrawCopy("COLZ");
 971
 972   canvas->cd(2);
 973   InitPadCOLZ();
 974   corrTrigger->DrawCopy("COLZ");
 975
 976   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
 977   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
 978
 979   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
 980   canvas->Divide(2, 1);
 981
 982   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
 983   corrVtx->GetYaxis()->SetRangeUser(0, 5);
 984
 985   canvas->cd(1);
 986   InitPadCOLZ();
 987   corrVtx->DrawCopy("COLZ");
 988
 989   canvas->cd(2);
 990   InitPadCOLZ();
 991   corrTrigger->DrawCopy("COLZ");
 992
 993   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
 994   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
 995 }
 996
 997 void TriggerBias(const char* fileName = "correction_map.root")
 998 {
 999   TFile* file = TFile::Open(fileName);
1000
1001   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
1002
1003   Prepare2DPlot(corr);
1004   corr->SetTitle("Trigger bias correction");
1005
1006   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
1007   InitPadCOLZ();
1008   corr->DrawCopy("COLZ");
1009
1010   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
1011   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
1012
1013   corr->GetYaxis()->SetRangeUser(0, 5);
1014
1015   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
1016   InitPadCOLZ();
1017   corr->DrawCopy("COLZ");
1018
1019   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
1020   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
1021 }
1022
1023 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1024 {
1025   gSystem->Load("libPWG0base");
1026
1027   TFile* file = TFile::Open(fileName);
1028   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1029   dNdEtaCorrection->LoadHistograms();
1030
1031   TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
1032   TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1033
1034   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
1035   canvas->Divide(2, 1);
1036
1037   canvas->cd(1);
1038   InitPad();
1039
1040   Prepare1DPlot(hist);
1041   hist->SetTitle("");
1042   hist->GetYaxis()->SetTitle("correction factor");
1043   hist->GetYaxis()->SetRangeUser(1, 1.5);
1044   hist->GetYaxis()->SetTitleOffset(1.6);
1045   hist->Draw();
1046
1047   canvas->cd(2);
1048   InitPad();
1049
1050   Prepare1DPlot(hist2);
1051   hist2->SetTitle("");
1052   hist2->GetYaxis()->SetTitle("correction factor");
1053   hist2->GetXaxis()->SetRangeUser(0, 5);
1054   hist2->GetYaxis()->SetTitleOffset(1.6);
1055   hist2->GetXaxis()->SetTitle("multiplicity");
1056   hist2->Draw();
1057
1058   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1059   pave->SetFillColor(0);
1060   pave->AddText("|z| < 10 cm");
1061   pave->Draw();
1062
1063   canvas->SaveAs("TriggerBias1D.eps");
1064 }
1065
1066 void VtxRecon()
1067 {
1068   TFile* file = TFile::Open("correction_map.root");
1069
1070   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
1071
1072   Prepare2DPlot(corr);
1073   corr->SetTitle("Vertex reconstruction correction");
1074
1075   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
1076   InitPadCOLZ();
1077   corr->DrawCopy("COLZ");
1078
1079   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
1080   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
1081
1082   corr->GetYaxis()->SetRangeUser(0, 5);
1083
1084   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
1085   InitPadCOLZ();
1086   corr->DrawCopy("COLZ");
1087
1088   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
1089   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
1090 }
1091
1092 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1093 {
1094   gSystem->Load("libPWG0base");
1095
1096   TFile* file = TFile::Open(fileName);
1097   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1098   dNdEtaCorrection->LoadHistograms();
1099
1100   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
1101   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1102
1103   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
1104   canvas->Divide(2, 1);
1105
1106   canvas->cd(1);
1107   InitPad();
1108
1109   Prepare1DPlot(hist);
1110   hist->SetTitle("");
1111   hist->GetYaxis()->SetTitle("correction factor");
1112   hist->GetYaxis()->SetRangeUser(1, 1.8);
1113   hist->GetYaxis()->SetTitleOffset(1.6);
1114   hist->DrawCopy();
1115
1116   canvas->cd(2);
1117   InitPad();
1118
1119   Prepare1DPlot(hist2);
1120   hist2->SetTitle("");
1121   hist2->GetYaxis()->SetTitle("correction factor");
1122   hist2->GetXaxis()->SetRangeUser(0, 20);
1123   hist2->GetYaxis()->SetTitleOffset(1.6);
1124   hist2->GetXaxis()->SetTitle("multiplicity");
1125   hist2->Draw();
1126
1127   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1128   pave->SetFillColor(0);
1129   pave->AddText("|z| < 10 cm");
1130   pave->Draw();
1131
1132   canvas->SaveAs("VtxRecon1D.eps");
1133
1134   Correction1DCreatePlots(fileName, folderName, 9.9, 2);
1135
1136   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1137   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1138
1139   Prepare1DPlot(corrX);
1140   Prepare1DPlot(corrZ);
1141
1142   corrX->GetYaxis()->SetTitleOffset(1.5);
1143   corrZ->GetYaxis()->SetTitleOffset(1.5);
1144
1145   corrX->SetTitle("a) z projection");
1146   corrZ->SetTitle("b) p_{T} projection");
1147
1148   corrX->GetYaxis()->SetTitle("Correction factor");
1149   corrZ->GetYaxis()->SetTitle("Correction factor");
1150
1151   corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
1152
1153   TString canvasName;
1154   canvasName.Form("VtxRecon1D_Track");
1155   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
1156   canvas->Divide(2, 1);
1157
1158   canvas->cd(1);
1159   InitPad();
1160   corrX->DrawCopy();
1161
1162   canvas->cd(2);
1163   InitPad();
1164   gPad->SetLogx();
1165   corrZ->Draw();
1166
1167   canvas->SaveAs("VtxRecon1D_Track.eps");
1168   canvas->SaveAs("VtxRecon1D_Track.gif");
1169 }
1170
1171 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
1172 {
1173   gSystem->Load("libPWG0base");
1174
1175   TFile::Open(fileName);
1176   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1177   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
1178
1179   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1180   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1181
1182   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1183   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1184
1185   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1186   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1187   gene->GetXaxis()->SetRangeUser(-10, 10);
1188   meas->GetXaxis()->SetRangeUser(-10, 10);
1189
1190   Float_t eff1 = gene->Integral() / meas->Integral();
1191   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1192
1193   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
1194
1195   gene->GetZaxis()->SetRangeUser(0.3, 10);
1196   meas->GetZaxis()->SetRangeUser(0.3, 10);
1197
1198   Float_t eff2 = gene->Integral() / meas->Integral();
1199   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1200
1201   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
1202
1203   gene->GetZaxis()->SetRangeUser(0.3, 1);
1204   meas->GetZaxis()->SetRangeUser(0.3, 1);
1205
1206   Float_t eff3 = gene->Integral() / meas->Integral();
1207   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1208
1209   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
1210 }
1211
1212 void Correction1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType = 0, Int_t correctionType2 = -1)
1213 {
1214   if (correctionType2 == -1)
1215     correctionType2 = correctionType;
1216
1217   TFile::Open(fileName);
1218   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1219   dNdEtaCorrection->LoadHistograms();
1220
1221   TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
1222   TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType2)->GetTrackCorrection()->GetMeasuredHistogram();
1223
1224   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1225   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1226   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1227   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1228   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kFALSE);
1229   gene->GetYaxis()->SetRange(0, 0);
1230   meas->GetYaxis()->SetRange(0, 0);
1231
1232   gene->GetXaxis()->SetRangeUser(-9.9, 9.9);
1233   meas->GetXaxis()->SetRangeUser(-9.9, 9.9);
1234   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kFALSE);
1235   gene->GetZaxis()->SetRange(0, 0);
1236   meas->GetZaxis()->SetRange(0, 0);
1237
1238   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1239   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1240   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kFALSE);
1241 }
1242
1243 TCanvas* Correction1D(Int_t correctionType = 0, const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType2 = -1)
1244 {
1245   gSystem->Load("libPWG0base");
1246
1247   Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType, correctionType2);
1248
1249   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1250   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1251   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1252
1253   Prepare1DPlot(corrX);
1254   Prepare1DPlot(corrY);
1255   Prepare1DPlot(corrZ);
1256
1257   /*
1258   corrX->SetTitle("a) z projection");
1259   corrY->SetTitle("b) #eta projection");
1260   corrZ->SetTitle("c) p_{T} projection");
1261   */
1262
1263   corrX->SetTitle("");
1264   corrY->SetTitle("");
1265   corrZ->SetTitle("");
1266
1267   corrX->SetTitleSize(0.06, "xyz");
1268   corrX->SetLabelSize(0.06, "xyz");
1269   corrY->SetTitleSize(0.06, "xyz");
1270   corrY->SetLabelSize(0.06, "xyz");
1271   corrZ->SetTitleSize(0.06, "xyz");
1272   corrZ->SetLabelSize(0.06, "xyz");
1273
1274   corrX->GetYaxis()->SetTitle("Correction factor");
1275   corrY->GetYaxis()->SetTitle("Correction factor");
1276   corrZ->GetYaxis()->SetTitle("Correction factor");
1277   //corrX->GetYaxis()->SetTitleOffset(1.7);
1278   //corrY->GetYaxis()->SetTitleOffset(1.7);
1279   //corrZ->GetYaxis()->SetTitleOffset(1.7);
1280   corrX->GetYaxis()->SetRangeUser(0.8, 1.5);
1281   corrY->GetYaxis()->SetRangeUser(0.8, 1.5);
1282   corrZ->GetYaxis()->SetRangeUser(0.8, 1.5);
1283
1284   corrZ->GetXaxis()->SetRangeUser(0.11, upperPtLimit);
1285
1286   TString canvasName;
1287   canvasName.Form(Form("Correction1D_%d_%s_%f", correctionType, fileName, upperPtLimit));
1288   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1289   canvas->Divide(3, 1);
1290
1291   TLatex* Tl = new TLatex;
1292   Tl->SetTextSize(0.06);
1293   Tl->SetBit(TLatex::kTextNDC);
1294
1295   canvas->cd(1);
1296   InitPad();
1297   gPad->SetTopMargin(0.05);
1298   gPad->SetBottomMargin(0.15);
1299   corrX->DrawCopy();
1300   Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1301   Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
1302
1303   canvas->cd(2);
1304   InitPad();
1305   gPad->SetTopMargin(0.05);
1306   gPad->SetBottomMargin(0.15);
1307   corrY->Draw();
1308   Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1309   Tl->DrawLatex(0.5, 0.8, "|vtx-z| < 10 cm");
1310
1311   canvas->cd(3);
1312   InitPad();
1313   gPad->SetTopMargin(0.05);
1314   gPad->SetBottomMargin(0.15);
1315   gPad->SetLogx();
1316   corrZ->Draw();
1317   corrZ->GetXaxis()->SetLabelOffset(0.005);
1318   corrZ->GetXaxis()->SetTitleOffset(1.2);
1319   Tl->DrawLatex(0.5, 0.88, "|vtx-z| < 10 cm");
1320   Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
1321
1322   return canvas;
1323 }
1324
1325 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1326 {
1327   gSystem->Load("libPWG0base");
1328
1329   Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
1330
1331   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1332   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1333   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1334
1335   Prepare1DPlot(corrX);
1336   Prepare1DPlot(corrY);
1337   Prepare1DPlot(corrZ);
1338
1339   corrX->SetTitle("a) z projection");
1340   corrY->SetTitle("a) #eta projection");
1341   corrZ->SetTitle("b) p_{T} projection");
1342
1343   corrY->GetYaxis()->SetTitle("correction factor");
1344   corrZ->GetYaxis()->SetTitle("correction factor");
1345
1346   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1347
1348   TString canvasName;
1349   canvasName.Form("Track2Particle1D_%s", folder);
1350   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1351   canvas->Divide(3, 1);
1352
1353   canvas->cd(1);
1354   InitPad();
1355   corrX->DrawCopy();
1356
1357   canvas->cd(2);
1358   InitPad();
1359   corrY->Draw();
1360
1361   canvas->cd(3);
1362   InitPad();
1363   corrZ->Draw();
1364
1365   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
1366   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
1367
1368   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
1369
1370   canvasName.Form("Track2Particle1D_%s_etapt", folder);
1371   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
1372   canvas->Divide(2, 1);
1373
1374   canvas->cd(1);
1375   InitPad();
1376   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
1377   corrY->GetYaxis()->SetRangeUser(1, 1.5);
1378   corrY->GetYaxis()->SetTitleOffset(1.5);
1379   corrY->DrawCopy();
1380   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
1381   pave->AddText("|z| < 10 cm");
1382   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
1383   pave->Draw();
1384
1385   canvas->cd(2);
1386   InitPad();
1387   gPad->SetLogx();
1388   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
1389   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
1390   corrZ->GetYaxis()->SetTitleOffset(1.5);
1391   corrZ->DrawCopy();
1392   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
1393   pave->AddText("|z| < 10 cm");
1394   pave->AddText("|#eta| < 0.8");
1395   pave->Draw();
1396
1397   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
1398   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
1399 }
1400
1401 /*
1402 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1403 {
1404   gSystem->Load("libPWG0base");
1405
1406   // particle type
1407   for (Int_t particle=0; particle<4; ++particle)
1408   {
1409     TString dirName;
1410     dirName.Form("correction_%d", particle);
1411     Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
1412
1413     TString tmpx, tmpy, tmpz;
1414     tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1415     tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1416     tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
1417
1418     TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1419     TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1420     TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
1421
1422     Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
1423
1424     TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1425     TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1426     TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
1427
1428     posX->Divide(negX);
1429     posY->Divide(negY);
1430     posZ->Divide(negZ);
1431
1432     Prepare1DPlot(posX);
1433     Prepare1DPlot(posY);
1434     Prepare1DPlot(posZ);
1435
1436     Float_t min = 0.8;
1437     Float_t max = 1.2;
1438
1439     posX->SetMinimum(min);
1440     posX->SetMaximum(max);
1441     posY->SetMinimum(min);
1442     posY->SetMaximum(max);
1443     posZ->SetMinimum(min);
1444     posZ->SetMaximum(max);
1445
1446     posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1447
1448     posX->GetYaxis()->SetTitleOffset(1.7);
1449     posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1450     posY->GetYaxis()->SetTitleOffset(1.7);
1451     posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1452     posZ->GetYaxis()->SetTitleOffset(1.7);
1453     posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
1454
1455     posZ->GetXaxis()->SetRangeUser(0, 1);
1456
1457     TString canvasName;
1458     canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
1459
1460     TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1461     canvas->Divide(3, 1);
1462
1463     canvas->cd(1);
1464     InitPad();
1465     posX->DrawCopy();
1466
1467     canvas->cd(2);
1468     InitPad();
1469     posY->DrawCopy();
1470
1471     canvas->cd(3);
1472     InitPad();
1473     posZ->DrawCopy();
1474
1475     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1476     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1477   }
1478 }
1479 */
1480
1481 void CompareTrack2Particle1D(const char* file1, const char* file2, Float_t upperPtLimit = 9.9)
1482 {
1483   loadlibs();
1484
1485   const char* folderName = "dndeta_correction";
1486
1487   c = new TCanvas("CompareTrack2Particle1D", "CompareTrack2Particle1D", 1200, 400);
1488   c->Divide(3, 1);
1489
1490   for (Int_t fileId = 0; fileId < 2; fileId++)
1491   {
1492     const char* file = ((fileId == 0) ? file1 : file2);
1493     Correction1DCreatePlots(file, folderName, upperPtLimit, 1);
1494
1495     TH1* corr[3];
1496     corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1497     corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y_div_measured_y"));
1498     corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1499     /*corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x"))->Clone(Form("hist_x_%d", fileId));
1500     corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y"))->Clone(Form("hist_y_%d", fileId));
1501     corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z"))->Clone(Form("hist_z_%d", fileId));*/
1502
1503     for (Int_t i=0; i<3; i++)
1504     {
1505       c->cd(i+1);
1506       InitPad();
1507       corr[i]->GetYaxis()->SetRangeUser(0.8, 2);
1508       corr[i]->SetLineColor(fileId+1);
1509       corr[i]->DrawCopy((fileId == 0) ? "" : "SAME");
1510     }
1511   }
1512
1513   return;
1514
1515   c->SaveAs(Form("%s.gif", canvas->GetName()));
1516   c->SaveAs(Form("%s.eps", canvas->GetName()));
1517 }
1518
1519 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
1520 {
1521   TFile::Open(fileName);
1522   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1523   dNdEtaCorrection->LoadHistograms();
1524
1525   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1526   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1527
1528   gene->GetZaxis()->SetRangeUser(0.2, 10);
1529   meas->GetZaxis()->SetRangeUser(0.2, 10);
1530   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
1531   gene->GetZaxis()->SetRange(0, 0);
1532   meas->GetZaxis()->SetRange(0, 0);
1533
1534   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1535   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1536   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
1537   gene->GetYaxis()->SetRange(0, 0);
1538   meas->GetYaxis()->SetRange(0, 0);
1539
1540   gene->GetXaxis()->SetRangeUser(-10, 10);
1541   meas->GetXaxis()->SetRangeUser(-10, 10);
1542   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
1543   gene->GetXaxis()->SetRange(0, 0);
1544   meas->GetXaxis()->SetRange(0, 0);
1545 }
1546
1547 TCanvas* Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1548 {
1549   gSystem->Load("libPWG0base");
1550
1551   Track2Particle2DCreatePlots(fileName);
1552
1553   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1554   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1555   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
1556
1557   Prepare2DPlot(corrYX);
1558   Prepare2DPlot(corrZX);
1559   Prepare2DPlot(corrZY);
1560
1561   const char* title = "";
1562   corrYX->SetTitle(title);
1563   corrZX->SetTitle(title);
1564   corrZY->SetTitle(title);
1565
1566   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1567   canvas->Divide(3, 1);
1568
1569   canvas->cd(1);
1570   InitPadCOLZ();
1571   corrYX->Draw("COLZ");
1572
1573   canvas->cd(2);
1574   InitPadCOLZ();
1575   corrZX->Draw("COLZ");
1576
1577   canvas->cd(3);
1578   InitPadCOLZ();
1579   corrZY->Draw("COLZ");
1580
1581   canvas->SaveAs(Form("corr_track2particle_%d.gif", gMax));
1582   canvas->SaveAs(Form("corr_track2particle_%d.eps", gMax));
1583
1584   return canvas;
1585 }
1586
1587 void CompareTrack2Particle2D()
1588 {
1589   gSystem->Load("libPWG0base");
1590
1591   Track2Particle2DCreatePlots("correction_maponly-positive.root");
1592
1593   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("pos_yx"));
1594   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("pos_zx"));
1595   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("pos_zy"));
1596
1597   Track2Particle2DCreatePlots("correction_maponly-negative.root");
1598
1599   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("neg_yx"));
1600   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("neg_zx"));
1601   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("neg_zy"));
1602
1603   posYX->Divide(negYX);
1604   posZX->Divide(negZX);
1605   posZY->Divide(negZY);
1606
1607   Prepare2DPlot(posYX);
1608   Prepare2DPlot(posZX);
1609   Prepare2DPlot(posZY);
1610
1611   Float_t min = 0.8;
1612   Float_t max = 1.2;
1613
1614   posYX->SetMinimum(min);
1615   posYX->SetMaximum(max);
1616   posZX->SetMinimum(min);
1617   posZX->SetMaximum(max);
1618   posZY->SetMinimum(min);
1619   posZY->SetMaximum(max);
1620
1621   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1622   canvas->Divide(3, 1);
1623
1624   canvas->cd(1);
1625   InitPadCOLZ();
1626   posYX->Draw("COLZ");
1627
1628   canvas->cd(2);
1629   InitPadCOLZ();
1630   posZX->Draw("COLZ");
1631
1632   canvas->cd(3);
1633   InitPadCOLZ();
1634   posZY->Draw("COLZ");
1635
1636   canvas->SaveAs("CompareTrack2Particle2D.gif");
1637   canvas->SaveAs("CompareTrack2Particle2D.eps");
1638 }
1639
1640 void Track2Particle3D()
1641 {
1642   // get left margin proper
1643
1644   TFile* file = TFile::Open("correction_map.root");
1645
1646   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1647
1648   corr->SetTitle("Correction Factor");
1649   SetRanges(corr->GetZaxis());
1650
1651   Prepare3DPlot(corr);
1652
1653   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1654   canvas->SetTheta(29.428);
1655   canvas->SetPhi(16.5726);
1656
1657   corr->Draw();
1658
1659   canvas->SaveAs("Track2Particle3D.gif");
1660   canvas->SaveAs("Track2Particle3D.eps");
1661 }
1662
1663 void Track2Particle3DAll()
1664 {
1665   TFile* file = TFile::Open("correction_map.root");
1666
1667   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1668   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1669   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1670
1671   gene->SetTitle("Generated Particles");
1672   meas->SetTitle("Measured Tracks");
1673   corr->SetTitle("Correction Factor");
1674
1675   Prepare3DPlot(gene);
1676   Prepare3DPlot(meas);
1677   Prepare3DPlot(corr);
1678
1679   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1680   canvas->Divide(3, 1);
1681
1682   canvas->cd(1);
1683   InitPad();
1684   gene->Draw();
1685
1686   canvas->cd(2);
1687   meas->Draw();
1688
1689   canvas->cd(3);
1690   corr->Draw();
1691
1692   canvas->SaveAs("Track2Particle3DAll.gif");
1693   canvas->SaveAs("Track2Particle3DAll.eps");
1694 }
1695
1696 void MultiplicityMC(Int_t xRangeMax = 50)
1697 {
1698   TFile* file = TFile::Open("multiplicityMC.root");
1699
1700   if (!file)
1701   {
1702     printf("multiplicityMC.root could not be opened.\n");
1703     return;
1704   }
1705
1706   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1707   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1708   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1709
1710   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1711   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1712   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1713   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1714   {
1715     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1716     proj->Fit("gaus", "0");
1717     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1718     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1719
1720     continue;
1721
1722     // draw for debugging
1723     new TCanvas;
1724     proj->DrawCopy();
1725     proj->GetFunction("gaus")->DrawCopy("SAME");
1726   }
1727
1728   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1729
1730   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1731   {
1732     Float_t mean = correction->GetBinContent(i);
1733     Float_t width = correctionWidth->GetBinContent(i);
1734
1735     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1736     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1737     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1738
1739     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1740     {
1741       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1742     }
1743   }
1744
1745   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1746   fMultiplicityESDCorrectedRebinned->Rebin(10);
1747   fMultiplicityESDCorrectedRebinned->Scale(0.1);
1748
1749   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1750   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1751   ratio->Divide(fMultiplicityMC);
1752
1753   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1754   ratio2->Divide(fMultiplicityMC);
1755
1756   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1757   canvas->Divide(3, 2);
1758
1759   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1760   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1761   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1762   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1763   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1764   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1765   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1766
1767   canvas->cd(1); //InitPad();
1768   fMultiplicityESD->Draw();
1769   fMultiplicityMC->SetLineColor(2);
1770   fMultiplicityMC->Draw("SAME");
1771
1772   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1773   legend->AddEntry(fMultiplicityESD, "ESD");
1774   legend->AddEntry(fMultiplicityMC, "MC");
1775   legend->Draw();
1776
1777   canvas->cd(2);
1778   fCorrelation->Draw("COLZ");
1779
1780   canvas->cd(3);
1781   correction->Draw();
1782   //correction->Fit("pol1");
1783   correctionWidth->SetLineColor(2);
1784   correctionWidth->Draw("SAME");
1785
1786   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1787   legend->AddEntry(correction, "#bar{x}");
1788   legend->AddEntry(correctionWidth, "#sigma");
1789   legend->Draw();
1790
1791   canvas->cd(4);
1792   ratio->Draw();
1793
1794   ratio2->SetLineColor(2);
1795   ratio2->Draw("SAME");
1796
1797   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1798   legend->AddEntry(ratio, "uncorrected");
1799   legend->AddEntry(ratio2, "corrected");
1800   legend->Draw();
1801
1802   canvas->cd(5);
1803   fMultiplicityESDCorrected->SetLineColor(kBlue);
1804   fMultiplicityESDCorrected->Draw();
1805   fMultiplicityMC->Draw("SAME");
1806   fMultiplicityESD->Draw("SAME");
1807
1808   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1809   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1810   legend->AddEntry(fMultiplicityMC, "MC");
1811   legend->AddEntry(fMultiplicityESD, "ESD");
1812   legend->Draw();
1813
1814   canvas->cd(6);
1815   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
1816   fMultiplicityESDCorrectedRebinned->Draw();
1817   fMultiplicityMC->Draw("SAME");
1818
1819   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1820   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1821   legend->AddEntry(fMultiplicityMC, "MC");
1822   legend->Draw();
1823
1824   canvas->SaveAs("MultiplicityMC.gif");
1825 }
1826
1827 void MultiplicityESD()
1828 {
1829   TFile* file = TFile::Open("multiplicityESD.root");
1830
1831   if (!file)
1832   {
1833     printf("multiplicityESD.root could not be opened.\n");
1834     return;
1835   }
1836
1837   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1838
1839   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1840
1841   fMultiplicityESD->Draw();
1842 }
1843
1844 void drawPlots(Int_t max)
1845 {
1846   gMax = max;
1847
1848   ptCutoff();
1849   TriggerBias();
1850   VtxRecon();
1851   Track2Particle2D();
1852   Track2Particle3D();
1853   ptSpectrum();
1854   dNdEta();
1855 }
1856
1857 void drawPlots()
1858 {
1859   drawPlots(5);
1860   drawPlots(2);
1861 }
1862
1863 void CompareCorrection2Measured(Float_t ptMin = 0.301, const char* dataInput = "analysis_esd_raw.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1864 {
1865   loadlibs();
1866
1867   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1868   TFile::Open(correctionMapFile);
1869   dNdEtaCorrection->LoadHistograms();
1870
1871   TFile* file = TFile::Open(dataInput);
1872
1873   if (!file)
1874   {
1875     cout << "Error. File not found" << endl;
1876     return;
1877   }
1878
1879   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1880   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1881
1882   gROOT->cd();
1883
1884   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
1885   hist1->SetTitle("mc");
1886   Printf("mc contains %f entries", hist1->Integral());
1887   Printf("mc contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), hist1->GetYaxis()->FindBin(-0.99), hist1->GetYaxis()->FindBin(0.99), hist1->GetZaxis()->FindBin(ptMin), hist1->GetNbinsZ()));
1888
1889   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
1890   hist2->SetTitle("esd");
1891   Printf("esd contains %f entries", hist2->Integral());
1892   Printf("esd contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), hist2->GetYaxis()->FindBin(-0.99), hist2->GetYaxis()->FindBin(0.99), hist2->GetZaxis()->FindBin(ptMin), hist2->GetNbinsZ()));
1893
1894   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1895   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
1896
1897   hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1898   hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1899   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
1900
1901   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
1902   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
1903   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
1904   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
1905   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
1906 }
1907
1908 void CompareCorrection2Generated(Float_t ptMin = 0.301, const char* dataInput = "analysis_mc.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1909 {
1910   loadlibs();
1911
1912   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1913   TFile::Open(correctionMapFile);
1914   dNdEtaCorrection->LoadHistograms();
1915
1916   TFile* file = TFile::Open(dataInput);
1917
1918   if (!file)
1919   {
1920     cout << "Error. File not found" << endl;
1921     return;
1922   }
1923
1924   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
1925   fdNdEtaAnalysis->LoadHistograms("dndetaTrVtx");
1926
1927   gROOT->cd();
1928
1929   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram()->Clone("mc");
1930   hist1->SetTitle("mc");
1931   Printf("mc contains %f entries", hist1->Integral());
1932   Printf("mc contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), hist1->GetYaxis()->FindBin(-0.99), hist1->GetYaxis()->FindBin(0.99), hist1->GetZaxis()->FindBin(ptMin), hist1->GetNbinsZ()));
1933
1934   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram()->Clone("esd");
1935   hist2->SetTitle("esd");
1936   Printf("esd contains %f entries", hist2->Integral());
1937   Printf("esd contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), hist2->GetYaxis()->FindBin(-0.99), hist2->GetYaxis()->FindBin(0.99), hist2->GetZaxis()->FindBin(ptMin), hist2->GetNbinsZ()));
1938
1939   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1940   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
1941
1942   hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1943   hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1944   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
1945
1946   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
1947   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
1948   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
1949   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
1950   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
1951 }
1952
1953 void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
1954 {
1955   loadlibs();
1956
1957   TFile* file = TFile::Open(dataInput);
1958
1959   if (!file)
1960   {
1961     cout << "Error. File not found" << endl;
1962     return;
1963   }
1964
1965   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1966   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1967
1968   TFile* file = TFile::Open(dataInput2);
1969
1970   if (!file)
1971   {
1972     cout << "Error. File not found" << endl;
1973     return;
1974   }
1975
1976   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
1977   fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
1978
1979   gROOT->cd();
1980
1981   TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
1982   hist1->SetTitle("esd1");
1983   Printf("esd1 contains %f entries", hist1->GetEntries());
1984   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()));
1985
1986   TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
1987   hist2->SetTitle("esd2");
1988   Printf("esd2 contains %f entries", hist2->GetEntries());
1989   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()));
1990
1991   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1992
1993   new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
1994   new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
1995   new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
1996
1997   TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
1998   TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
1999
2000   Printf("event1 contains %f entries", event1->GetEntries());
2001   Printf("event2 contains %f entries", event2->GetEntries());
2002   Printf("event1 integral is %f", event1->Integral());
2003   Printf("event2 integral is %f", event2->Integral());
2004   Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
2005   Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
2006
2007   projx1 = event1->ProjectionX();
2008   projx2 = event2->ProjectionX();
2009
2010   new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
2011
2012   projx1->Divide(projx2);
2013   new TCanvas; projx1->Draw();
2014
2015   event1->Divide(event2);
2016   new TCanvas; event1->Draw("COLZ");
2017
2018 }
2019
2020 void DrawTrackletOrigin(const char* fileName = "correction_map.root", Bool_t myFile = kTRUE)
2021 {
2022   TFile::Open(fileName);
2023
2024   Int_t maxHists = 8;
2025   TH1* hist[8];
2026
2027   const Int_t kRebin = 8;
2028
2029   const char* titles[] = { "PP", "SS", "PP'", "PS'", "PS", "SP'", "SS'", "" };
2030
2031   if (myFile)
2032   {
2033     for (Int_t i=0; i<maxHists; i++)
2034     {
2035       hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
2036       if (hist[i]->GetDimension() == 2)
2037         hist[i] = ((TH2*) hist[i])->ProjectionX(Form("fDeltaPhi_clone_%d", i));
2038     }
2039   }
2040   else
2041   {
2042     maxHists = 6;
2043     const char* names[] = { "DePhiPPTracklets", "DePhiSecTracklets", "DePhiPpTracklets", "DePhiPSTracklets", "DePhiPSdaugTracklets", "DePhiSPTracklets" };
2044     for (Int_t i=0; i<maxHists; i++)
2045       hist[i] = (TH1*) gFile->Get(names[i]);
2046   }
2047
2048   // clone before rebinning
2049   good = (TH1*) hist[0]->Clone("good");
2050   good->Add(hist[4]);
2051
2052   bad = (TH1*) hist[1]->Clone("bad");
2053   bad->Add(hist[2]);
2054   bad->Add(hist[3]);
2055   bad->Add(hist[5]);
2056   if (myFile)
2057     bad->Add(hist[6]);
2058
2059   c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c");
2060   TH1* ref = 0;
2061   Bool_t nw = kFALSE;
2062   if (!c)
2063   {
2064     c = new TCanvas("c", "c", 600, 600);
2065     nw = kTRUE;
2066     ref = (TH1*) c->GetListOfPrimitives()->At(1);
2067   }
2068   c->cd();
2069   c->SetRightMargin(0.05);
2070   c->SetTopMargin(0.05);
2071   c->SetLogy();
2072   c->SetGridx();
2073   c->SetGridy();
2074
2075   Int_t order[] = { 0, 4, 1, 2, 3, 5, 6, 7 };
2076   //Int_t colors[]  = {1,2,4,1,2,4,1,2,4};
2077   Int_t colors[]  = {1,2,3,4,6,7,8,102};
2078   Int_t markers[]  = {20, 21, 22, 23, 24, 25, 26, 27, 28};
2079
2080   TLegend* legend = new TLegend(0.75, 0.6, 0.93, 0.93);
2081   legend->SetFillColor(0);
2082   legend->SetTextSize(0.04);
2083
2084   Int_t total = 0;
2085   for (Int_t ii=0; ii<maxHists; ii++)
2086   {
2087     i = order[ii];
2088
2089     hist[i]->Rebin(kRebin);
2090     hist[i]->SetStats(kFALSE);
2091     hist[i]->SetLineColor(colors[i]);
2092     hist[i]->SetLineWidth(2);
2093     //hist[i]->SetMarkerStyle(markers[i]);
2094     //hist[i]->SetMarkerColor(colors[i]);
2095     //hist[i]->SetLineStyle(ii+1);
2096     hist[i]->GetXaxis()->SetRangeUser(-0.09, 0.09);
2097     hist[i]->GetYaxis()->SetRangeUser(5, hist[i]->GetMaximum() * 2);
2098     hist[i]->GetYaxis()->SetTitleOffset(1.3);
2099     hist[i]->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2100
2101     if (i == 0 && ref)
2102       hist[i]->Scale(1.0 / hist[i]->GetMaximum() * ref->GetMaximum());
2103
2104     hist[i]->DrawCopy(((i == 0 && nw) ? "" : "SAME"));
2105
2106     total += hist[i]->GetEntries();
2107
2108     if (i != 7)
2109       legend->AddEntry(hist[i], titles[i], "L");
2110   }
2111
2112   legend->Draw();
2113   c->SaveAs("spd_tracklets_deltaphi_detailed.eps");
2114
2115   Printf("Total: %d", total);
2116   for (Int_t i=0; i<maxHists; i++)
2117     Printf("Histogram %d (%s) contains %.2f %% of the entries", i, titles[i], 100.0 * hist[i]->GetEntries() / total);
2118
2119   printf("|  Delta phi  |  Acc. %%  |  ");
2120   for (Int_t i=0; i<maxHists; i++)
2121     printf("%3s %%   |  ", titles[i]);
2122   Printf("");
2123
2124   for (Float_t f = 0.01; f < 0.09; f += 0.01)
2125   {
2126     Int_t integralBegin = hist[0]->GetXaxis()->FindBin(-f);
2127     Int_t integralEnd = hist[0]->GetXaxis()->FindBin(f);
2128
2129     Int_t total2 = 0;
2130     for (Int_t i=0; i<maxHists; i++)
2131       total2 += (Int_t) hist[i]->Integral(integralBegin, integralEnd);
2132
2133     printf("|    %.2f     |  %6.2f  |  ", f, 100.0 * total2 / total);
2134
2135     for (Int_t i=0; i<maxHists; i++)
2136       printf("%6.2f  |  ", (hist[i]->GetEntries() > 0) ? (100.0 * hist[i]->Integral(integralBegin, integralEnd) / hist[i]->GetEntries()) : -1.0);
2137     Printf("");
2138   }
2139
2140   eff = new TH1F("eff", ";#Delta#varphi cut (rad.)", 101,-0.0005, 0.1005);
2141   cont = new TH1F("cont", "cont", 101,-0.0005, 0.1005);
2142   signalOverBg = new TH1F("signalOverBg", "signalOverBg", 101,-0.0005, 0.1005);
2143   for (Float_t cut=0.000; cut<0.10; cut += 0.001)
2144   {
2145     Float_t accGood = good->Integral(good->GetXaxis()->FindBin(-cut), good->GetXaxis()->FindBin(cut));
2146     Float_t accBad = bad->Integral(bad->GetXaxis()->FindBin(-cut), bad->GetXaxis()->FindBin(cut));
2147     Float_t sB = accGood / accBad;
2148     eff->Fill(cut, 100.0 * accGood / good->Integral());
2149     cont->Fill(cut, 100.0 * accBad / (accGood + accBad));
2150     signalOverBg->Fill(cut, sB);
2151   }
2152
2153   //new TCanvas; signalOverBg->Draw();
2154
2155   c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c2");
2156   Bool_t nw = kFALSE;
2157   if (!c)
2158   {
2159     c = new TCanvas("c2", "c2", 600, 600);
2160     nw = kTRUE;
2161   }
2162   c->cd();
2163   c->SetRightMargin(0.05);
2164   c->SetTopMargin(0.05);
2165   c->SetGridx();
2166   c->SetGridy();
2167   gPad->SetLogy();
2168   good->Rebin(kRebin);
2169   bad->Rebin(kRebin);
2170   good->GetXaxis()->SetRangeUser(-0.09, 0.09);
2171   good->GetYaxis()->SetTitleOffset(1.3);
2172   good->SetStats(0);
2173   good->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2174   good->DrawCopy((nw) ? "" : "SAME");
2175
2176   bad->SetLineColor(2);
2177   bad->SetLineStyle(2);
2178   bad->SetLineWidth(2);
2179   //bad->SetMarkerColor(2);
2180   //bad->SetMarkerStyle(7);
2181   bad->DrawCopy("SAME");
2182
2183   TLegend* legend = new TLegend(0.2, 0.13, 0.85, 0.25);
2184   legend->SetFillColor(0);
2185   legend->SetTextSize(0.04);
2186   legend->AddEntry(good, "Primaries", "L");
2187   legend->AddEntry(bad, "Secondaries + Background", "L");
2188   legend->Draw();
2189
2190   c->SaveAs("spd_tracklets_deltaphi.eps");
2191
2192   c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c3");
2193   Bool_t nw = kFALSE;
2194   if (!c)
2195   {
2196     c = new TCanvas("c3", "c3", 600, 600);
2197     nw = kTRUE;
2198   }
2199   c->cd();
2200   c->SetRightMargin(0.05);
2201   c->SetTopMargin(0.05);
2202   c->SetGridx();
2203   c->SetGridy();
2204
2205   TLegend* legend = new TLegend(0.5, 0.6, 0.93, 0.75);
2206   legend->SetFillColor(0);
2207   legend->SetTextSize(0.04);
2208   legend->AddEntry(eff, "Efficiency (%)", "L");
2209   legend->AddEntry(cont, "Contamination (%)", "L");
2210
2211   eff->SetStats(0);
2212   eff->GetXaxis()->SetRangeUser(0, 0.08);
2213   eff->GetYaxis()->SetRangeUser(1e-3, 105);
2214   eff->SetLineWidth(2);
2215   eff->DrawCopy((nw) ? "" : "SAME");
2216   cont->SetLineStyle(2);
2217   cont->SetLineWidth(2);
2218   cont->SetLineColor(2);
2219   cont->DrawCopy("SAME");
2220   legend->Draw();
2221
2222   c->SaveAs("spd_tracklets_efficiency.eps");
2223
2224
2225 }
2226
2227 void Tracklets_Asymmetry()
2228 {
2229   TFile::Open("correction_map.root");
2230
2231   Int_t maxHists = 7;
2232   TH1* hist[8];
2233
2234   Int_t colors[]  = {1,2,3,4,6,7,8,102};
2235   const char* titles[] = { "PP", "SS", "PP'", "PS'", "PS", "SP'", "SS'", "" };
2236
2237   TLegend* legend = new TLegend(0.75, 0.6, 0.93, 0.93);
2238
2239   for (Int_t i=0; i<maxHists; i++)
2240   {
2241     hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
2242     hist[i]->Rebin(10);
2243
2244     for (Int_t j=hist[i]->GetNbinsX()/2; j<=hist[i]->GetNbinsX(); j++)
2245       if (hist[i]->GetBinContent(j) > 0)
2246         hist[i]->SetBinContent(j, (hist[i]->GetBinContent(j) -  hist[i]->GetBinContent(hist[i]->GetXaxis()->FindBin(-hist[i]->GetXaxis()->GetBinCenter(j)))) / hist[i]->GetBinContent(j));
2247
2248     hist[i]->SetStats(kFALSE);
2249     hist[i]->SetLineColor(colors[i]);
2250     hist[i]->GetXaxis()->SetRangeUser(0.001, 0.09);
2251     //hist[i]->GetYaxis()->SetRangeUser(5, hist[i]->GetMaximum() * 2);
2252     hist[i]->GetYaxis()->SetTitleOffset(1.3);
2253     hist[i]->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2254     hist[i]->Draw(((i == 0) ? "" : "SAME"));
2255
2256     legend->AddEntry(hist[i], titles[i], "L");
2257   }
2258
2259   legend->Draw();
2260 }
2261
2262 TH2* GetCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2263 {
2264   // returns the correction factor with pt integrated out
2265
2266   loadlibs();
2267
2268   TFile::Open(fileName);
2269
2270   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
2271   if (!dNdEtaCorrection->LoadHistograms())
2272     return;
2273
2274   //  hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
2275
2276   gener = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
2277   measu = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
2278
2279   gener->GetZaxis()->SetRange(gener->GetZaxis()->FindBin(ptmin), gener->GetNbinsZ()+1);
2280   TH2D *gener_xy = gener->Project3D("yx");
2281
2282   measu->GetZaxis()->SetRange(measu->GetZaxis()->FindBin(ptmin), measu->GetNbinsZ()+1);
2283   TH2D *measu_xy = measu->Project3D("yx");
2284
2285   cout << measu->GetZaxis()->FindBin(ptmin) << " " << measu->GetNbinsZ()+1 << endl;
2286
2287   TCanvas *canp = new TCanvas("canp","canp",600,1000);
2288   canp->Divide(1,2,0.0001,0.0001);
2289   canp->cd(1);
2290   gener_xy->Draw("COLZ");
2291   canp->cd(2);
2292   measu_xy->Draw("COLZ");
2293
2294
2295   TCanvas *canpr = new TCanvas("canpr","canpr",700,500);
2296   canpr->cd();
2297   TH2D *proj = new TH2D(*gener_xy);
2298   proj->Divide(measu_xy);
2299
2300 //   proj = hist->Project3D("yx");
2301   proj->Draw("COLZ");
2302
2303   return proj;
2304 }
2305
2306 void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2307 {
2308   TH2* proj = GetCorrection(fileName, dirName, ptmin);
2309
2310   const Float_t limit = 5;
2311
2312   TString array = "{";
2313   TString arrayEnd = "}";
2314
2315   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2316   {
2317     Int_t begin = -1;
2318     Int_t end = -1;
2319     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2320     {
2321       if (begin == -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2322         begin = x;
2323       if (begin != -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2324         end = x;
2325     }
2326     Printf("Limits for y = %d are %d to %d", y, begin, end);
2327
2328     if (y > 1)
2329       array += ", ";
2330     array += Form("%d", begin);
2331
2332     if (y > 1)
2333       arrayEnd.Prepend(", ");
2334     arrayEnd.Prepend(Form("%d", (end == -1) ? -1 : proj->GetNbinsX() + 1 - end));
2335   }
2336   array += "}";
2337   arrayEnd.Prepend("{");
2338
2339   Printf("Begin array:");
2340   Printf("%s", array.Data());
2341
2342   Printf("End array (mirrored) (should be the same):");
2343   Printf("%s", arrayEnd.Data());
2344 }
2345
2346 void AverageMultiplicity(const char* fileName = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
2347 {
2348   loadlibs();
2349
2350   TFile::Open(fileName);
2351
2352   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
2353   dNdEtaCorrection->LoadHistograms();
2354   TH2* events = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram();
2355   TH3* tracks = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram();
2356
2357   Float_t nEvents = events->Integral(events->GetXaxis()->FindBin(-1), events->GetXaxis()->FindBin(1), 0, events->GetNbinsY()+1);
2358   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);
2359
2360   Printf("%f %f --> %f", nEvents, nTracks, nTracks / nEvents);
2361 }
2362
2363 void GetAverageCorrectionFactor(Float_t etaRange = 1.5, Float_t vertexRange = 9.9, const char* rawFile = "analysis_esd_raw.root", const char* mcFile = "analysis_mc.root")
2364 {
2365   loadlibs();
2366
2367   TFile::Open(rawFile);
2368   dNdEtaAnalysis* raw = new dNdEtaAnalysis("dndeta", "dndeta");
2369   raw->LoadHistograms("fdNdEtaAnalysisESD");
2370   raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->GetXaxis()->SetRangeUser(-vertexRange, vertexRange);
2371   tracks = raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Project3D("y");
2372   events = raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("events", 0, raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->GetNbinsY() + 1);
2373   Float_t nEvents = events->Integral(events->FindBin(-vertexRange), events->FindBin(vertexRange));
2374   tracks->Scale(1.0 / nEvents / tracks->GetBinWidth(1));
2375
2376   TFile::Open(mcFile);
2377   dNdEtaAnalysis* mc = new dNdEtaAnalysis("dndeta", "dndeta");
2378   mc->LoadHistograms("dndetaTrVtx");
2379   mcH = mc->GetdNdEtaPtCutOffCorrectedHistogram(0);
2380
2381   new TCanvas;
2382   mcH->SetLineColor(2);
2383   mcH->DrawCopy();
2384   tracks->DrawCopy("SAME");
2385
2386   new TCanvas;
2387   mcH->GetYaxis()->SetRangeUser(0, 5);
2388   mcH->Divide(tracks);
2389   mcH->DrawCopy();
2390   mcH->Fit("pol0", "", "", -etaRange, etaRange);
2391 }
2392
2393 void TrackCuts_Comparison(char* histName, Int_t plotWhich = 0, const char* fileName = "correction_map.root", Bool_t mirror = kFALSE)
2394 {
2395   // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
2396   //    --> manually disable it in the run.C
2397   //
2398   // plotWhich: 0 = only before
2399   //            1 = both
2400   //            2 = only after
2401   //
2402   // mirror: kTRUE --> project negative values on the positive side
2403
2404
2405   file = TFile::Open(fileName);
2406
2407   Int_t count = 0;
2408   Int_t colors[] = { 1, 2, 3, 4, 5, 6 };
2409
2410   TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
2411   TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
2412   TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
2413
2414   TCanvas* c1 = new TCanvas(histName, histName, 800, 1200);
2415   c1->Divide(1, 2);
2416   //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
2417   //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
2418
2419   const char* folders2[] = { "before_cuts", "after_cuts" };
2420   Bool_t first = kTRUE;
2421   for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
2422   {
2423     const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
2424     const char* names[] =    { "all", "primaries", "secondaries" };
2425     TH1* base = 0;
2426     TH1* prim = 0;
2427     TH1* sec = 0;
2428     for (Int_t i = 0; i < 3; i++)
2429     {
2430       TString folder;
2431       folder.Form("%s/%s/%s", folders1[i], folders2[j], histName);
2432       TH1* hist = (TH1*) file->Get(folder);
2433
2434       if (mirror)
2435       {
2436         for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
2437         {
2438           Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
2439           if (bin != newBin)
2440           {
2441             hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
2442             hist->SetBinContent(bin, 0);
2443           }
2444         }
2445       }
2446
2447       legend->AddEntry(hist, Form("%s %s", names[i], folders2[j]));
2448
2449       c1->cd(1);
2450       hist->SetLineColor(colors[count]);
2451       hist->DrawCopy((count == 0) ? "" : "SAME");
2452
2453       switch (i)
2454       {
2455         case 0: base = hist; break;
2456         case 1: prim = hist; break;
2457         case 2: sec = hist; break;
2458       }
2459
2460       count++;
2461     }
2462
2463     TH1* eff    = (TH1*) prim->Clone("eff"); eff->Reset();
2464     TH1* purity = (TH1*) prim->Clone("purity"); purity->Reset();
2465
2466     for (Int_t bin = 1; bin <= prim->GetNbinsX(); bin++)
2467     {
2468       eff->SetBinContent(bin, prim->Integral(1, bin) / prim->Integral(1, prim->GetNbinsX() + 1));
2469       if (prim->Integral(1, bin) + sec->Integral(1, bin) > 0)
2470         purity->SetBinContent(bin, sec->Integral(1, bin) / (prim->Integral(1, bin) + sec->Integral(1, bin)));
2471     }
2472
2473     eff->GetYaxis()->SetRangeUser(0, 1);
2474     eff->SetLineColor(colors[0+j*2]);
2475     eff->SetStats(kFALSE);
2476     purity->SetLineColor(colors[1+j*2]);
2477
2478     legend3->AddEntry(eff, Form("%s: efficiency", folders2[j]));
2479     legend3->AddEntry(purity, Form("%s: contamination", folders2[j]));
2480
2481     c1->cd(2);
2482     eff->DrawCopy((first) ? "" : "SAME");
2483     first = kFALSE;
2484     purity->DrawCopy("SAME");
2485   }
2486
2487   c1->cd(1)->SetLogy();
2488   c1->cd(1)->SetGridx();
2489   c1->cd(1)->SetGridy();
2490   legend->Draw();
2491
2492   //c2->cd();
2493  // c2->SetGridx();
2494  // c2->SetGridy();
2495   //legend2->Draw();
2496
2497   c1->cd(2)->SetGridx();
2498   c1->cd(2)->SetGridy();
2499   legend3->Draw();
2500
2501   //c1->SaveAs(Form("%s.png", histName));
2502 }
2503
2504 void TrackCuts_DCA()
2505 {
2506   file = TFile::Open("correction_map.root");
2507   hist = (TH2*) file->Get("esd_track_cuts/before_cuts/dXYvsDZ");
2508
2509   TCanvas* c1 = new TCanvas("c1", "c1", 600, 600);
2510   c1->SetLogz();
2511   c1->SetRightMargin(0.12);
2512   c1->SetBottomMargin(0.12);
2513
2514   hist->SetStats(kFALSE);
2515   hist->Draw("COLZ");
2516
2517   ellipse = new TEllipse(0, 0, 4);
2518   ellipse->SetLineWidth(2);
2519   ellipse->SetLineStyle(2);
2520   ellipse->SetFillStyle(0);
2521   ellipse->Draw();
2522
2523   c1->SaveAs("trackcuts_dca_2d.eps");
2524 }
2525
2526 void FindNSigma(TH2* hist, Int_t nSigma = 1)
2527 {
2528   TH1* proj = hist->ProjectionY();
2529   proj->Reset();
2530
2531   for (Int_t bin=1; bin<=proj->GetNbinsX(); bin++)
2532   {
2533     if (hist->Integral(1, hist->GetNbinsX(), bin, bin) == 0)
2534       continue;
2535
2536     Int_t limit = -1;
2537     for (limit = 1; limit<=hist->GetNbinsX(); limit++)
2538   }
2539 }
2540
2541 void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2542 {
2543   TH2* proj = GetCorrection(fileName, dirName, ptmin);
2544
2545   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2546     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2547       if (proj->GetBinContent(x, y) > 5 || proj->GetBinContent(x, y) == 0)
2548       {
2549         proj->SetBinContent(x, y, 0);
2550       }
2551       else
2552         proj->SetBinContent(x, y, 1);
2553
2554
2555   input->Multiply(proj);
2556 }
2557
2558 void MakeGaussianProfile(const char* histName = "fVertexCorrelation", Bool_t subtractMean = kFALSE)
2559 {
2560     TFile::Open("correction_map.root");
2561
2562     TH2* hist2d = (TH2*) gFile->Get(histName);
2563     hist2d->Sumw2();
2564
2565     TH1* result = hist2d->ProjectionX("result");
2566     result->GetYaxis()->SetTitle(hist2d->GetYaxis()->GetTitle());
2567     result->Reset();
2568
2569     for (Int_t x=1; x<hist2d->GetNbinsX(); ++x)
2570     {
2571         hist = hist2d->ProjectionY(Form("temp_%d", x), x, x);
2572         if (hist->GetEntries() == 0)
2573             continue;
2574         if (hist->Fit("gaus") == 0)
2575         {
2576             func = hist->GetFunction("gaus");
2577             mean = func->GetParameter(1);
2578             error = func->GetParError(1);
2579
2580             if (subtractMean)
2581                 mean = hist2d->GetXaxis()->GetBinCenter(x) - mean;
2582
2583             result->SetBinContent(x, mean);
2584             result->SetBinError(x, error);
2585
2586             if (x % 10 == 0)
2587             {
2588                 new TCanvas;
2589                 ((TH1*) hist->Clone())->DrawCopy();
2590             }
2591         }
2592         //break;
2593     }
2594
2595     new TCanvas;
2596     result->GetYaxis()->SetRangeUser(-0.2, 0.2);
2597     result->Draw();
2598 }
2599
2600 TH2* GetAcceptance(void* corr2d_void)
2601 {
2602         corr2d = (AliCorrectionMatrix2D*) corr2d_void;
2603         corr_xy = (TH2*) corr2d->GetCorrectionHistogram()->Clone("acceptance");
2604
2605         // fold in acceptance
2606         for (Int_t x=1; x<=corr_xy->GetNbinsX(); ++x)
2607                 for (Int_t y=1; y<=corr_xy->GetNbinsY(); ++y)
2608                 {
2609                         if (corr_xy->GetBinContent(x, y) > 1.5)
2610                                 corr_xy->SetBinContent(x, y, 0);
2611
2612                         if (corr_xy->GetBinContent(x, y) > 0)
2613                                 corr_xy->SetBinContent(x, y, 1);
2614
2615                         corr_xy->SetBinError(x, y, 0);
2616                 }
2617
2618         return corr_xy;
2619 }
2620
2621 void ZeroOutsideAcceptance(TH2* acc, TH3* hist)
2622 {
2623   for (Int_t x=0; x<=acc->GetNbinsX()+1; ++x)
2624     for (Int_t y=0; y<=acc->GetNbinsY()+1; ++y)
2625     {
2626       if (acc->GetBinContent(x, y) > 1.5 || acc->GetBinContent(x, y) == 0)
2627       {
2628         for (Int_t z=0; z<=hist->GetNbinsZ()+1; ++z)
2629         {
2630           hist->SetBinContent(x, y, z, 0);
2631           hist->SetBinError(x, y, z, 0);
2632         }
2633       }
2634     }
2635 }
2636
2637 void DrawPhi()
2638 {
2639   loadlibs();
2640
2641   TFile::Open("correction_map.root");
2642   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
2643   if (!dNdEtaCorrection->LoadHistograms())
2644     return 0;
2645
2646   TFile::Open("analysis_esd.root");
2647   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
2648   fdNdEtaAnalysis->LoadHistograms();
2649
2650   // acc. map!
2651   //acc = GetAcceptance(dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000));
2652   acc = dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000)->GetCorrectionHistogram();
2653   //new TCanvas; acc->Draw("COLZ");
2654
2655   histG = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
2656   ZeroOutsideAcceptance(acc, histG);
2657   //new TCanvas; histG->Project3D("yx")->Draw("COLZ");
2658   //histG->GetYaxis()->SetRangeUser(-0.9, 0.9);
2659
2660   histM = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
2661   ZeroOutsideAcceptance(acc, histM);
2662   //histM->GetYaxis()->SetRangeUser(-0.9, 0.9);
2663
2664   TFile::Open("analysis_mc.root");
2665   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndetaTrVtxMC", "dndetaTrVtxMC");
2666   fdNdEtaAnalysis2->LoadHistograms("dndetaTrVtx");
2667
2668   histMC = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
2669   ZeroOutsideAcceptance(acc, histMC);
2670   //new TCanvas; histMC->Project3D("yx2")->Draw("COLZ");
2671
2672   //histG->GetZaxis()->SetRangeUser(1,2); histMC->GetZaxis()->SetRangeUser(1,2);
2673   new TCanvas; a = histG->Project3D("yx3"); a->Add(histMC->Project3D("yx4"), -1); a->Draw("COLZ");
2674
2675   //histMC->GetYaxis()->SetRangeUser(-0.9, 0.9);
2676
2677   c = new TCanvas;
2678
2679   histG->GetXaxis()->SetRangeUser(-9.9, 9.9);
2680   histG->Project3D("z")->Draw();
2681
2682   histM->GetXaxis()->SetRangeUser(-9.9, 9.9);
2683   proj = histM->Project3D("z2");
2684   proj->SetLineColor(2);
2685   proj->Draw("SAME");
2686
2687   histMC->GetXaxis()->SetRangeUser(-9.9, 9.9);
2688   projMC = histMC->Project3D("z3");
2689   projMC->SetLineColor(3);
2690   projMC->Draw("SAME");
2691 }
2692
2693 void PrintEventStats(Int_t corrID = 3)
2694 {
2695   loadlibs();
2696
2697   /*
2698   TFile::Open("analysis_mc.root");
2699   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD");
2700   fdNdEtaAnalysis->LoadHistograms();
2701   trackHist = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
2702   eventHist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetGeneratedHistogram();
2703   */
2704
2705   TFile::Open("correction_map.root");
2706   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
2707   if (!dNdEtaCorrection->LoadHistograms())
2708     return;
2709   trackHist = dNdEtaCorrection->GetCorrection(corrID)->GetTrackCorrection()->GetGeneratedHistogram();
2710   eventHist = dNdEtaCorrection->GetCorrection(corrID)->GetEventCorrection()->GetGeneratedHistogram();
2711
2712   trackHist->GetXaxis()->SetRange(0, trackHist->GetNbinsX()+1);
2713   trackHist->GetZaxis()->SetRange(0, trackHist->GetNbinsZ()+1);
2714   eta = trackHist->Project3D("y");
2715
2716   events = eventHist->Integral(0, eventHist->GetNbinsX()+1, 0, eventHist->GetNbinsY()+1);
2717
2718   eta->Scale(1.0 / events);
2719
2720   Float_t avgN = eta->Integral(0, eta->GetNbinsX()+1);
2721   Printf("<N> = %f", avgN);
2722
2723   eta->Scale(1.0 / eta->GetXaxis()->GetBinWidth(1));
2724
2725   Printf("dndeta | eta = 0 is %f", (eta->GetBinContent(eta->FindBin(0.01)) + eta->GetBinContent(eta->FindBin(-0.01))) / 2);
2726   Printf("dndeta in |eta| < 1 is %f", eta->Integral(eta->FindBin(-0.99), eta->FindBin(0.99)) / (eta->FindBin(0.99) - eta->FindBin(-0.99) + 1));
2727   Printf("dndeta in |eta| < 1.5 is %f", eta->Integral(eta->FindBin(-1.49), eta->FindBin(1.49)) / (eta->FindBin(1.49) - eta->FindBin(-1.49) + 1));
2728
2729   stats = (TH2*) gFile->Get("fEventStats");
2730   proj = stats->ProjectionX();
2731   gROOT->ProcessLine(".L PrintHist.C");
2732   PrintHist2D(stats);
2733   PrintHist(proj);
2734
2735   Printf("+++ TRIGGER EFFICIENCIES +++");
2736
2737   Printf("INEL = %.1f", 100. * (proj->GetBinContent(1) - stats->GetBinContent(1, 1)) / proj->GetBinContent(1));
2738   Printf("NSD  = %.1f", 100. * (proj->GetBinContent(2) - stats->GetBinContent(2, 1)) / proj->GetBinContent(2));
2739   Printf("ND  = %.1f",  100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1)) / proj->GetBinContent(3));
2740   Printf("SD  = %.1f",  100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1)) / proj->GetBinContent(4));
2741   Printf("DD  = %.1f",  100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1)) / proj->GetBinContent(5));
2742
2743   for (Int_t i=7; i<=proj->GetNbinsX(); i++)
2744     if (proj->GetBinContent(i) > 0)
2745       Printf("bin %d (process type %d) = %.2f", i, (Int_t) proj->GetXaxis()->GetBinCenter(i), 100.0 * (proj->GetBinContent(i) - stats->GetBinContent(i, 1)) / proj->GetBinContent(i));
2746
2747   //eta->Draw();
2748 }
2749
2750 void TestAsymmetry()
2751 {
2752   loadlibs();
2753
2754   TFile* file2 = TFile::Open("analysis_mc.root");
2755
2756   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
2757   fdNdEtaAnalysis->LoadHistograms();
2758   fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
2759
2760   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy();
2761
2762   hist = (TH1*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
2763   hist2 = (TH1*) hist->Clone("hist2");
2764   for (Int_t x=1; x<=hist->GetNbinsX(); x++)
2765     for (Int_t y=1; y<=hist->GetNbinsY(); y++)
2766       for (Int_t z=1; z<=hist->GetNbinsZ(); z++)
2767       {
2768         Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
2769         hist->SetBinContent(x, y, z, hist2->GetBinContent(hist->GetNbinsX() / 2, TMath::Min(y, hist->GetNbinsY() + 1 - y), z));
2770       }
2771
2772   hist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
2773   for (Int_t x=1; x<=hist->GetNbinsX(); x++)
2774     for (Int_t y=1; y<=hist->GetNbinsY(); y++)
2775       {
2776         //Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
2777         hist->SetBinContent(x, y, hist->GetBinContent(hist->GetNbinsX() / 2, y));
2778       }
2779
2780   fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
2781   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerColor(2);
2782   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetLineColor(2);
2783   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerStyle(5);
2784   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy("SAMEP");
2785 }
2786
2787 void DeltaPhiFromPt(Float_t smearing = 0.005)
2788 {
2789   loadlibs();
2790
2791   TFile::Open("analysis_mc.root");
2792   hist = (TH1*) gFile->Get("dndeta_check_pt");
2793
2794   dPhiHist = new TH1F("dPhiHist", ";#Delta phi", 400, -0.1, 0.1);
2795   dPhiHist2 = new TH1F("dPhiHist2", ";#Delta phi", 400, -0.1, 0.1);
2796
2797   for (Int_t i=1; i<=hist->GetNbinsX(); i++)
2798   {
2799     Float_t pt = hist->GetBinCenter(i);
2800     Float_t deltaPhi = (0.076 - 0.039) / 2 / (pt / 0.15);
2801
2802     if (smearing > 0)
2803     {
2804       gaus = new TF1("mygaus", "gaus(0)", -0.1, 0.1);
2805       gaus->SetParameters(1, -deltaPhi, smearing);
2806
2807       dPhiHist->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
2808
2809       dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
2810       gaus->SetParameters(1, deltaPhi, smearing);
2811       dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
2812     }
2813     else
2814 {
2815 dPhiHist->Fill(deltaPhi, hist->GetBinContent(i) / 2);
2816 dPhiHist2->Fill(deltaPhi, hist->GetBinContent(i) / 2);
2817 dPhiHist2->Fill(-deltaPhi, hist->GetBinContent(i) / 2);
2818 }
2819   }
2820
2821   new TCanvas;
2822   dPhiHist->Draw();
2823   dPhiHist2->SetLineColor(2);
2824   dPhiHist2->Draw("SAME");
2825   gPad->SetLogy();
2826
2827   TFile::Open("trackletsDePhi.root");
2828   //TFile::Open("tmp/correction_maponly-positive.root");
2829   //TFile::Open("tmp/correction_map.root");
2830   //tracklets = (TH1*) gFile->Get(Form("fDeltaPhi_%d", 0));
2831   tracklets = (TH1*) gFile->Get("DePhiPPTracklets");
2832   tracklets->Scale(1.0 / tracklets->GetMaximum() * dPhiHist->GetMaximum());
2833   tracklets->SetLineColor(4);
2834   tracklets->Draw("SAME");
2835 }