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 || 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.2);
 149     return;
 150   }
 151 }
 152
 153 void PrintAllInfos()
 154 {
 155   PrintInfo();
 156
 157   Printf("RAW ESD");
 158   TFile::Open("analysis_esd_raw.root");
 159   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 160   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
 161   fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
 162
 163   const Int_t num = 3;
 164   const char* results[] = { "dndeta", "dndetaTr", "dndetaTrVtx" };
 165
 166   TFile::Open("analysis_esd.root");
 167   for (Int_t i=0; i<num; i++)
 168   {
 169     Printf("CORRECTED %s", results[i]);
 170     dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 171     fdNdEtaAnalysis->LoadHistograms(results[i]);
 172     fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
 173   }
 174 }
 175
 176 void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
 177 {
 178   gSystem->Load("libPWG0base");
 179
 180   TFile::Open(esdFile);
 181   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 182   fdNdEtaAnalysisESD->LoadHistograms();
 183
 184   TFile::Open(mcFile);
 185   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
 186   fdNdEtaAnalysisMC->LoadHistograms();
 187   //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 188
 189   for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
 190     fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
 191
 192   fdNdEtaAnalysisESD->DrawHistograms();
 193 }
 194
 195 void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 196 {
 197   gSystem->Load("libPWG0base");
 198
 199   const char* ESDfolder = 0;
 200
 201   if (plot == 0) // all
 202     ESDfolder = "dndeta";
 203   else if (plot == 1) // mb
 204     ESDfolder = "dndeta_mb";
 205   else if (plot == 2) // mb vtx
 206     ESDfolder = "dndeta_mbvtx";
 207
 208   TFile::Open("analysis_esd.root");
 209   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 210   fdNdEtaAnalysisESD->LoadHistograms();
 211
 212   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 213   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 214
 215   TH2F* hist = 0;
 216
 217   if (plot == 0) // all
 218     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
 219   else if (plot == 1) // mb
 220     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
 221   else if (plot == 2) // mb vtx
 222     hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
 223
 224   TH1* proj = hist->ProjectionX();
 225
 226   TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
 227   for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
 228   {
 229     Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
 230     if (value != 0)
 231     {
 232       printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
 233       vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
 234     }
 235   }
 236
 237   new TCanvas;
 238   vertex->DrawCopy();
 239 }
 240
 241 void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 242 {
 243   loadlibs();
 244
 245   TFile::Open("analysis_esd.root");
 246   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
 247   fdNdEtaAnalysisESD->LoadHistograms();
 248
 249   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 250   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 251
 252   //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
 253   //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 254
 255   TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
 256   TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 257
 258   TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
 259
 260   new TCanvas;
 261   histESD->Draw();
 262
 263   new TCanvas;
 264   histCorr->Draw();
 265
 266   for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
 267     for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
 268       for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
 269       {
 270         Float_t value1 = histESD->GetBinContent(x, y, z);
 271         Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
 272
 273         if (value2 > 0 && value1 > 0)
 274         {
 275           printf("%f %f %f\n", value1, value2, value1 / value2);
 276           diff->Fill(value1 / value2);
 277         }
 278       }
 279
 280   new TCanvas;
 281   diff->Draw();
 282 }
 283
 284 Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
 285 {
 286   Double_t avgMC = 0;
 287   Double_t avgESD = 0;
 288   for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
 289   {
 290     avgMC += histMC->GetBinContent(bin);
 291     avgESD += histESD->GetBinContent(bin);
 292   }
 293   Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
 294
 295   avgMC /= nBins;
 296   avgESD /= nBins;
 297
 298   // deviation when integrate in |eta| < 0.8 between mc and esd
 299   Double_t diffFullRange = (avgMC - avgESD) / avgMC;
 300
 301   Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
 302
 303   return diffFullRange;
 304 }
 305
 306 void dNdEtaNoResolution()
 307 {
 308   loadlibs();
 309
 310   TFile::Open("correction_map.root");
 311
 312   const char* correctionMapFolder = "dndeta_correction";
 313   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 314   dNdEtaCorrection->LoadHistograms();
 315   dNdEtaCorrection->GetTrack2ParticleCorrection()->PrintInfo(0.3);
 316
 317   TFile::Open("analysis_mc.root");
 318   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
 319   fdNdEtaAnalysis->LoadHistograms();
 320   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, "ESD (no resolution effect) -> MB with vertex");
 321   fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0)->SetMarkerStyle(21);
 322
 323   TFile::Open("analysis_mc.root");
 324   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
 325   fdNdEtaAnalysisMC->LoadHistograms();
 326   fdNdEtaAnalysisMC->Finish(0, 0, AlidNdEtaCorrection::kNone, "MC: MB with vertex");
 327
 328   DrawdNdEtaRatio(fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0), fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(0), "MB with vertex (no resolution effect)", 3);
 329 }
 330
 331 TH1* GetMCHist(const char* folder, Float_t ptCut, const char* tag)
 332 {
 333   loadlibs();
 334
 335   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(folder, folder);
 336   fdNdEtaAnalysis->LoadHistograms();
 337   fdNdEtaAnalysis->Finish(0, ptCut, AlidNdEtaCorrection::kNone, tag);
 338   return fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 339 }
 340
 341 void dNdEtaFinal(Bool_t spd = kTRUE)
 342 {
 343   TFile* file = TFile::Open("analysis_esd.root");
 344   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
 345   TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
 346   Prepare1DPlot(histESD);
 347   Prepare1DPlot(histESDnsd);
 348
 349   TCanvas* canvas = new TCanvas("dNdEtaFinal", "dNdEtaFinal", 600, 600);
 350   gPad->SetTopMargin(0.05);
 351   gPad->SetRightMargin(0.05);
 352   gPad->SetLeftMargin(0.12);
 353   gPad->SetBottomMargin(0.12);
 354   gPad->SetGridx();
 355   gPad->SetGridy();
 356
 357   Float_t etaMax = 1.9;
 358   Float_t histMax = 1.39;
 359   Float_t systErrorValue = 0.023;
 360   Float_t systErrorNSDValue = 0.081;
 361   if (!spd)
 362   {
 363     //etaMax = 1.5;
 364     histMax = 0.99;
 365     systErrorValue = 0.043;
 366     systErrorNSDValue = 0.088;
 367   }
 368
 369   dummy = new TH2F("dummy", ";#eta;dN_{ch}/d#eta", 100, -etaMax, etaMax, 100, 3, 8);
 370   dummy->SetStats(0);
 371   dummy->GetYaxis()->SetTitleOffset(1.3);
 372
 373   histESD->SetMarkerStyle(20);
 374   histESDnsd->SetMarkerStyle(21);
 375   histESDnsd->SetMarkerColor(4);
 376   histESDnsd->SetLineColor(4);
 377   histESD->SetMarkerSize(1.5);
 378   histESDnsd->SetMarkerSize(1.5);
 379
 380   histESD->GetXaxis()->SetRangeUser(-histMax, histMax);
 381   histESDnsd->GetXaxis()->SetRangeUser(-histMax, histMax);
 382
 383   legend = new TLegend(0.3, 0.2, 0.78, 0.4);
 384   legend->SetFillColor(0);
 385   legend->SetTextSize(0.04);
 386   legend->AddEntry(histESD, "Inelastic events", "P");
 387   legend->AddEntry(histESDnsd, "NSD events", "P");
 388
 389   dummy->Draw();
 390
 391   // syst errors.
 392   TH1* systError = (TH1*) histESD->Clone("systError");
 393   for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
 394     systError->SetBinError(i, systError->GetBinContent(i) * systErrorValue);
 395   // change error drawing style
 396   systError->SetFillColor(15);
 397   systError->DrawCopy("SAME E2 ][");
 398
 399   // syst error NSD
 400   for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
 401   {
 402     systError->SetBinContent(i, histESDnsd->GetBinContent(i));
 403     systError->SetBinError(i, systError->GetBinContent(i) * systErrorNSDValue);
 404   }
 405   systError->DrawCopy("SAME E2 ][");
 406
 407   histESD->Draw("SAME");
 408   histESDnsd->Draw("SAME");
 409   legend->Draw();
 410
 411   canvas->SaveAs(Form("%s_dndeta_final.eps", (spd) ? "spd" : "tpc"));
 412 }
 413
 414 void dNdEtaPythiaPhojet()
 415 {
 416   // evtl. deactivate acceptance maps in dNdEtaAnalysis.cxx
 417
 418   loadlibs();
 419
 420   TH1* hist[4];
 421
 422   TFile::Open("LHC08c11_10TeV_0.5T/mb1/spd/analysis_mc.root");
 423   hist[0] =         (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
 424   hist[1] =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
 425
 426   TFile::Open("LHC08c15_10TeV_0.5T_Phojet/mb1/spd/analysis_mc.root");
 427   hist[2] =         (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMCPhojet");
 428   hist[3] =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsdPhojet");
 429
 430   file = TFile::Open("pythia_phojet_dndeta.root", "RECREATE");
 431   for (Int_t i=0; i<4; i++)
 432     hist[i]->Write();
 433   file->Close();
 434 }
 435
 436 void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 437 {
 438   loadlibs();
 439
 440   TFile* file = TFile::Open("analysis_esd.root");
 441
 442   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 443   fdNdEtaAnalysis->LoadHistograms("dndeta");
 444
 445   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
 446   TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
 447   TH1* histESDnsdNoPt = (TH1*) file->Get("dndetaNSD/dNdEta");
 448   TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
 449   TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
 450   TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
 451   TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
 452   TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
 453   TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
 454
 455   Prepare1DPlot(histESD);
 456   Prepare1DPlot(histESDnsd);
 457   Prepare1DPlot(histESDMB);
 458   Prepare1DPlot(histESDMBVtx);
 459
 460   Prepare1DPlot(histESDNoPt);
 461   Prepare1DPlot(histESDMBNoPt);
 462   Prepare1DPlot(histESDMBVtxNoPt);
 463   Prepare1DPlot(histESDMBTracksNoPt);
 464
 465   histESD->SetLineWidth(0);
 466   histESDnsd->SetLineWidth(0);
 467   histESDMB->SetLineWidth(0);
 468   histESDMBVtx->SetLineWidth(0);
 469
 470   histESDNoPt->SetLineWidth(0);
 471   histESDMBNoPt->SetLineWidth(0);
 472   histESDMBVtxNoPt->SetLineWidth(0);
 473
 474   histESD->SetMarkerColor(1);
 475   histESDnsd->SetMarkerColor(6);
 476   histESDMB->SetMarkerColor(2);
 477   histESDMBVtx->SetMarkerColor(4);
 478
 479   histESD->SetLineColor(1);
 480   histESDnsd->SetLineColor(6);
 481   histESDMB->SetLineColor(2);
 482   histESDMBVtx->SetLineColor(4);
 483
 484   histESDNoPt->SetMarkerColor(1);
 485   histESDMBNoPt->SetMarkerColor(2);
 486   histESDMBVtxNoPt->SetMarkerColor(4);
 487   histESDMBTracksNoPt->SetMarkerColor(3);
 488
 489   histESD->SetMarkerStyle(20);
 490   histESDnsd->SetMarkerStyle(29);
 491   histESDMB->SetMarkerStyle(21);
 492   histESDMBVtx->SetMarkerStyle(22);
 493
 494   histESDNoPt->SetMarkerStyle(20);
 495   histESDMBNoPt->SetMarkerStyle(21);
 496   histESDMBVtxNoPt->SetMarkerStyle(22);
 497   histESDMBTracksNoPt->SetMarkerStyle(23);
 498
 499   Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 0.89 : 1.79;
 500   Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 1.2 : 2.3;
 501   //Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 0.89 : 1.39;
 502   //Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 1.2 : 1.9;
 503
 504   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 505   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 506   histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 507   histESDnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 508
 509   histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 510   histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 511   histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 512   histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 513
 514   Float_t max = TMath::Max(histESDMBVtx->GetMaximum(), histESDMB->GetMaximum());
 515   max = TMath::Max(max, histESD->GetMaximum());
 516   max = TMath::Max(max, histESDMBTracksNoPt->GetMaximum());
 517
 518   TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.4);
 519   legend->SetFillColor(0);
 520   legend->AddEntry(histESDMBVtx, "Triggered, vertex");
 521   legend->AddEntry(histESDMB, "Triggered");
 522   legend->AddEntry(histESD, "All events");
 523
 524   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, max * 1.1);
 525   dummy->GetYaxis()->SetRangeUser(2.1, max * 1.1);
 526   Prepare1DPlot(dummy);
 527   dummy->SetStats(kFALSE);
 528   dummy->SetXTitle("#eta");
 529   dummy->SetYTitle("dN_{ch}/d#eta");
 530   dummy->GetYaxis()->SetTitleOffset(1);
 531
 532   TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
 533
 534   dummy->DrawCopy();
 535   histESDMBVtx->Draw("SAME");
 536   histESDMB->Draw("SAME");
 537   histESD->Draw("SAME");
 538   legend->Draw();
 539
 540   if (save)
 541   {
 542     canvas->SaveAs("dNdEta1.gif");
 543     canvas->SaveAs("dNdEta1.eps");
 544   }
 545
 546   if (onlyESD)
 547     return;
 548
 549   loadlibs();
 550
 551   TFile* file2 = TFile::Open("analysis_mc.root");
 552
 553   TH1* histMCTrVtx =       (TH1*) GetMCHist("dndetaTrVtx", -1, "MC: MB with trigger")->Clone("histMCTrVtx");
 554   TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
 555
 556   TH1* histMC =            (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
 557   TH1* histMCTr =          (TH1*) GetMCHist("dndetaTr", -1, "MC: minimum bias")->Clone("histMCTr");
 558   TH1* histMCnsd =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
 559
 560   TH1* histMCPtCut =       (TH1*) GetMCHist("dndeta", 0.21, "MC: full inelastic, pt cut")->Clone("histMCPtCut");
 561   TH1* histMCTrPtCut =     (TH1*) GetMCHist("dndetaTr", 0.21, "MC: minimum bias, pt cut")->Clone("histMCTrPtCut");
 562   TH1* histMCTrVtxPtCut =  (TH1*) GetMCHist("dndetaTrVtx", 0.21, "MC: MB with trigger, pt cut")->Clone("histMCTrVtxPtCut");
 563   TH1* histMCnsdNoPt =     (TH1*) GetMCHist("dndetaNSD", 0.21, "MC: NSD, put cut")->Clone("histMCnsdNoPt");
 564   TH1* histMCTracksPtCut = (TH1*) GetMCHist("dndetaTracks", 0.21, "MC: Tracks w/o resolution effect, pt cut")->Clone("histMCTracksPtCut");
 565
 566   Prepare1DPlot(histMC);
 567   Prepare1DPlot(histMCnsd);
 568   Prepare1DPlot(histMCTr);
 569   Prepare1DPlot(histMCTrVtx);
 570
 571   Prepare1DPlot(histMCPtCut);
 572   Prepare1DPlot(histMCTrPtCut);
 573   Prepare1DPlot(histMCTrVtxPtCut);
 574   Prepare1DPlot(histMCTracksPtCut);
 575
 576   histMC->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 577   histMCnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 578   histMCTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 579   histMCTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 580
 581   histMCPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 582   histMCTrPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 583   histMCTrVtxPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 584   histMCTracksPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 585
 586   histMC->SetLineColor(1);
 587   histMCnsd->SetLineColor(6);
 588   histMCTr->SetLineColor(2);
 589   histMCTrVtx->SetLineColor(4);
 590
 591   histMCPtCut->SetLineColor(1);
 592   histMCTrPtCut->SetLineColor(2);
 593   histMCTrVtxPtCut->SetLineColor(4);
 594   if (histMCTracksPtCut)
 595     histMCTracksPtCut->SetLineColor(3);
 596
 597   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
 598
 599   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
 600   dummy2->GetYaxis()->SetRangeUser(0, max * 1.1);
 601
 602   dummy2->DrawCopy();
 603   histMC->Draw("SAME");
 604   histMCnsd->Draw("SAME");
 605   histMCTr->Draw("SAME");
 606   histMCTrVtx->Draw("SAME");
 607   histESD->Draw("SAME");
 608   histESDnsd->Draw("SAME");
 609   histESDMB->Draw("SAME");
 610   histESDMBVtx->Draw("SAME");
 611   histESDNoPt->Draw("SAME");
 612   histESDMBNoPt->Draw("SAME");
 613   histESDMBVtxNoPt->Draw("SAME");
 614   histESDMBTracksNoPt->Draw("SAME");
 615   histMCPtCut->Draw("SAME");
 616   histMCTrPtCut->Draw("SAME");
 617   histMCTrVtxPtCut->Draw("SAME");
 618   if (histMCTracksPtCut)
 619     histMCTracksPtCut->Draw("SAME");
 620
 621   if (save)
 622   {
 623     canvas2->SaveAs("dNdEta2.gif");
 624     canvas2->SaveAs("dNdEta2.eps");
 625   }
 626
 627   TH1* ratio = (TH1*) DrawdNdEtaRatio(histESD, histMC, "full_inelastic", etaPlotLimit)->Clone();
 628   TH1* ratioTr = (TH1*) DrawdNdEtaRatio(histESDMB, histMCTr, "triggered", etaPlotLimit)->Clone();
 629   TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
 630   TH1* ratioTrVtxNoPt = (TH1*) DrawdNdEtaRatio(histESDMBVtxNoPt, histMCTrVtxPtCut, "triggered_vertex_nopt", etaPlotLimit)->Clone();
 631   TH1* ratioNSD = (TH1*) DrawdNdEtaRatio(histESDnsd, histMCnsd, "NSD", etaPlotLimit)->Clone();
 632
 633   // draw ratios of single steps
 634   c7 = new TCanvas("all_ratios", "all_ratios", 600, 600);
 635   c7->SetRightMargin(0.05);
 636   c7->SetTopMargin(0.05);
 637   c7->SetGridx();
 638   c7->SetGridy();
 639
 640   ratioTrVtxNoPt->SetMarkerStyle(20);
 641   ratioTrVtx->SetMarkerStyle(21);
 642   ratioTr->SetMarkerStyle(23);
 643   ratio->SetMarkerStyle(22);
 644   ratioNSD->SetMarkerStyle(26);
 645
 646   ratioTrVtxNoPt->SetMarkerSize(2);
 647   ratioTrVtx->SetMarkerSize(2);
 648   ratioTr->SetMarkerSize(2);
 649   ratio->SetMarkerSize(2);
 650   ratioNSD->SetMarkerSize(2);
 651
 652   ratioTrVtxNoPt->SetMarkerColor(1);
 653   ratioTrVtx->SetMarkerColor(2);
 654   ratioTr->SetMarkerColor(4);
 655   ratio->SetMarkerColor(2);
 656   ratioNSD->SetMarkerColor(1);
 657
 658   ratioTrVtxNoPt->SetLineColor(1);
 659   ratioTrVtx->SetLineColor(2);
 660   ratioTr->SetLineColor(4);
 661   ratio->SetLineColor(2);
 662   ratioNSD->SetLineColor(1);
 663
 664   legend7 = new TLegend(0.13, 0.7, 0.94, 0.9);
 665   legend7->SetFillColor(0);
 666   legend7->SetTextSize(0.035);
 667   legend7->SetNColumns(2);
 668
 669   flat = new TF1("flat", "-1", -5, 5);
 670   ratioTrVtxNoPt->Add(flat);
 671   ratioTrVtx->Add(flat);
 672   ratioTr->Add(flat);
 673   ratio->Add(flat);
 674   ratioNSD->Add(flat);
 675
 676   ratioTrVtxNoPt->Scale(100);
 677   ratioTrVtx->Scale(100);
 678   ratioTr->Scale(100);
 679   ratio->Scale(100);
 680   ratioNSD->Scale(100);
 681
 682   ratio->Add(ratioTr, -1);
 683   ratioNSD->Add(ratioTr, -1);
 684   ratioTr->Add(ratioTrVtx, -1);
 685   ratioTrVtx->Add(ratioTrVtxNoPt, -1);
 686
 687   legend7->AddEntry(ratioTrVtxNoPt, "Track-to-particle", "P");
 688   legend7->AddEntry(ratio, "Trigger-bias INEL", "P");
 689   legend7->AddEntry(ratioTr, "Vertex-reconstruction", "P");
 690   legend7->AddEntry(ratioNSD, "Trigger-bias NSD", "P");
 691   if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
 692     legend7->AddEntry(ratioTrVtx, "p_{T} cut-off", "P");
 693
 694   TH1* dummy7 = new TH2F("dummy7", ";#eta;Deviation in %", 100, -etaPlotLimit, etaPlotLimit, 100, -5, 7);
 695   dummy7->SetStats(0);
 696   dummy7->Draw();
 697
 698   ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 699   ratioTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 700   ratioTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 701   ratioTrVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 702   ratioNSD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 703
 704   ratio->Draw("HIST EP SAME");
 705   ratioTr->Draw("HIST EP SAME");
 706   if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
 707     ratioTrVtx->Draw("HIST EP SAME");
 708   ratioTrVtxNoPt->Draw("HIST EP SAME");
 709   ratioNSD->Draw("HIST EP SAME");
 710   legend7->Draw();
 711
 712   c7->SaveAs("ratios.eps");
 713
 714   new TCanvas;
 715   dummy2->DrawCopy();
 716   histMCnsd->Draw("SAME");
 717   histESDnsd->Draw("SAME");
 718
 719   ratio = (TH1*) histMC->Clone("ratio");
 720   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
 721
 722   ratio->Divide(histESD);
 723   ratioNoPt->Divide(histESDNoPt);
 724
 725   ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 726
 727   ratio->SetLineColor(1);
 728   ratioNoPt->SetLineColor(2);
 729
 730   Double_t average = 0;       // average deviation from 1 in ratio (depends on the number of bins if statistical)
 731   for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
 732     average += TMath::Abs(ratio->GetBinContent(bin) - 1);
 733   Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
 734   average /= nBins;
 735   Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
 736
 737   PrintIntegratedDeviation(histMC, histESD, "all events");
 738   PrintIntegratedDeviation(histMCnsd, histESDnsd, "all events (NSD)");
 739   PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
 740   PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
 741   PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
 742   PrintIntegratedDeviation(histMCnsdNoPt, histESDnsdNoPt, "all events (NSD) (no pt corr)");
 743   PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
 744   PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
 745   PrintIntegratedDeviation(histESD, histESDNoPt, "pt cut off correction");
 746
 747   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 600, 600);
 748   canvas3->Range(0, 0, 1, 1);
 749   //canvas3->Divide(1, 2, 0, 0);
 750
 751   //canvas3->cd(1);
 752   TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
 753   pad1->SetTopMargin(0.05);
 754   pad1->SetLeftMargin(0.13);
 755   pad1->Draw();
 756
 757   TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
 758   pad2->SetLeftMargin(0.13);
 759   pad2->Draw();
 760
 761   pad1->SetRightMargin(0.01);
 762   pad2->SetRightMargin(0.01);
 763
 764   // no border between them
 765   pad1->SetBottomMargin(0);
 766   pad2->SetTopMargin(0);
 767
 768   pad1->cd();
 769   pad1->SetGridx();
 770   pad1->SetGridy();
 771
 772   legend->AddEntry(histMC, "MC prediction");
 773
 774   dummy->GetXaxis()->SetLabelSize(0.08);
 775   dummy->GetYaxis()->SetLabelSize(0.08);
 776   dummy->GetXaxis()->SetTitleSize(0.08);
 777   dummy->GetYaxis()->SetTitleSize(0.08);
 778   dummy->GetYaxis()->SetTitleOffset(0.8);
 779   dummy->DrawCopy();
 780   histESDMBVtx->Draw("SAME");
 781   histESDMB->Draw("SAME");
 782   histESD->Draw("SAME");
 783   histMC->Draw("SAME");
 784
 785   legend->SetTextSize(0.08);
 786   legend->Draw();
 787
 788   pad2->cd();
 789   pad2->SetBottomMargin(0.15);
 790   //pad2->SetGridx();
 791   //pad2->SetGridy();
 792
 793   Float_t minR = 0.91; //TMath::Min(0.961, ratio->GetMinimum() * 0.95);
 794   Float_t maxR = 1.09; //TMath::Max(1.049, ratio->GetMaximum() * 1.05);
 795
 796   TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
 797   dummy3.SetStats(kFALSE);
 798   for (Int_t i=1; i<=100; ++i)
 799     dummy3.SetBinContent(i, 1);
 800   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
 801   dummy3.SetLineWidth(2);
 802   dummy3.GetXaxis()->SetLabelSize(0.08);
 803   dummy3.GetYaxis()->SetLabelSize(0.08);
 804   dummy3.GetXaxis()->SetTitleSize(0.08);
 805   dummy3.GetYaxis()->SetTitleSize(0.08);
 806   dummy3.GetYaxis()->SetTitleOffset(0.8);
 807   dummy3.DrawCopy();
 808
 809   ratio->Draw("SAME");
 810
 811   //pad2->Draw();
 812
 813   canvas3->Modified();
 814
 815   if (save)
 816   {
 817     canvas3->SaveAs("dNdEta.gif");
 818     canvas3->SaveAs("dNdEta.eps");
 819   }
 820
 821   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
 822
 823   ratio->Draw();
 824   ratioNoPt->Draw("SAME");
 825
 826   TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
 827   legend->SetFillColor(0);
 828   legend->AddEntry(ratio, "mc/esd");
 829   legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
 830   legend->Draw();
 831 }
 832
 833 TH1* DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
 834 {
 835   TCanvas* canvas3 = new TCanvas(name, name, 600, 600);
 836   canvas3->Range(0, 0, 1, 1);
 837
 838   TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
 839   pad1->Draw();
 840
 841   TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
 842   pad2->Draw();
 843
 844   pad1->SetRightMargin(0.01);
 845   pad2->SetRightMargin(0.01);
 846   pad1->SetTopMargin(0.05);
 847   pad1->SetLeftMargin(0.13);
 848   pad2->SetLeftMargin(0.13);
 849   pad2->SetBottomMargin(0.15);
 850
 851   // no border between them
 852   pad1->SetBottomMargin(0);
 853   pad2->SetTopMargin(0);
 854
 855   pad1->cd();
 856   pad1->SetGridx();
 857   pad1->SetGridy();
 858
 859   TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.3);
 860   legend->SetFillColor(0);
 861   legend->AddEntry(corr, "Corrected");
 862   legend->AddEntry(mc, "MC prediction");
 863   legend->SetTextSize(0.08);
 864
 865   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 2.7, corr->GetMaximum() * 1.1);
 866   Prepare1DPlot(dummy);
 867   dummy->SetStats(kFALSE);
 868   dummy->SetXTitle("#eta");
 869   dummy->SetYTitle("dN_{ch}/d#eta");
 870   dummy->GetYaxis()->SetTitleOffset(1);
 871
 872   dummy->GetXaxis()->SetLabelSize(0.08);
 873   dummy->GetYaxis()->SetLabelSize(0.08);
 874   dummy->GetXaxis()->SetTitleSize(0.08);
 875   dummy->GetYaxis()->SetTitleSize(0.08);
 876   dummy->GetYaxis()->SetTitleOffset(0.8);
 877   dummy->DrawCopy();
 878
 879   corr->Draw("SAME");
 880   mc->Draw("SAME");
 881
 882   legend->Draw();
 883
 884   pad2->cd();
 885   pad2->SetBottomMargin(0.15);
 886   //pad2->SetGridx();
 887   //pad2->SetGridy();
 888
 889   TH1* ratio = (TH1*) mc->Clone("ratio");
 890   ratio->Divide(corr);
 891
 892   Float_t minR = TMath::Min(0.91, ratio->GetMinimum() * 0.95);
 893   Float_t maxR = TMath::Max(1.09, ratio->GetMaximum() * 1.05);
 894
 895   TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
 896   dummy3.SetStats(kFALSE);
 897   for (Int_t i=1; i<=100; ++i)
 898         dummy3.SetBinContent(i, 1);
 899   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
 900   dummy3.SetLineWidth(2);
 901   dummy3.GetXaxis()->SetLabelSize(0.08);
 902   dummy3.GetYaxis()->SetLabelSize(0.08);
 903   dummy3.GetXaxis()->SetTitleSize(0.08);
 904   dummy3.GetYaxis()->SetTitleSize(0.08);
 905   dummy3.GetYaxis()->SetTitleOffset(0.8);
 906   dummy3.DrawCopy();
 907
 908   ratio->Draw("SAME");
 909
 910   canvas3->Modified();
 911
 912   return ratio;
 913 }
 914
 915 void ptSpectrum()
 916 {
 917   TFile* file = TFile::Open("analysis_esd.root");
 918   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
 919
 920   TFile* file2 = TFile::Open("analysis_mc.root");
 921   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
 922
 923   TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
 924   InitPad();
 925   gPad->SetLogy();
 926
 927   Prepare1DPlot(histMC);
 928   Prepare1DPlot(histESD);
 929
 930   histESD->SetTitle("");
 931   histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 932   histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
 933
 934   histMC->SetLineColor(kBlue);
 935   histESD->SetLineColor(kRed);
 936
 937   histESD->GetYaxis()->SetTitleOffset(1.5);
 938   histESD->GetXaxis()->SetRangeUser(0, 4.9999);
 939
 940   histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
 941
 942   histESD->Draw();
 943   histMC->Draw("SAME");
 944
 945   canvas->SaveAs("ptSpectrum.gif");
 946   canvas->SaveAs("ptSpectrum.eps");
 947 }
 948
 949 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
 950 {
 951   gSystem->Load("libPWG0base");
 952
 953   TFile::Open(fileName);
 954   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
 955   dNdEtaCorrection->LoadHistograms();
 956
 957   TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
 958   TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
 959
 960   Prepare2DPlot(corrTrigger);
 961   corrTrigger->SetTitle("b) Trigger bias correction");
 962
 963   Prepare2DPlot(corrVtx);
 964   corrVtx->SetTitle("a) Vertex reconstruction correction");
 965
 966   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
 967   corrVtx->GetYaxis()->SetTitle("Multiplicity");
 968
 969   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
 970   canvas->Divide(2, 1);
 971
 972   canvas->cd(1);
 973   InitPadCOLZ();
 974   corrVtx->DrawCopy("COLZ");
 975
 976   canvas->cd(2);
 977   InitPadCOLZ();
 978   corrTrigger->DrawCopy("COLZ");
 979
 980   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
 981   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
 982
 983   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
 984   canvas->Divide(2, 1);
 985
 986   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
 987   corrVtx->GetYaxis()->SetRangeUser(0, 5);
 988
 989   canvas->cd(1);
 990   InitPadCOLZ();
 991   corrVtx->DrawCopy("COLZ");
 992
 993   canvas->cd(2);
 994   InitPadCOLZ();
 995   corrTrigger->DrawCopy("COLZ");
 996
 997   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
 998   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
 999 }
1000
1001 void TriggerBias(const char* fileName = "correction_map.root")
1002 {
1003   TFile* file = TFile::Open(fileName);
1004
1005   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
1006
1007   Prepare2DPlot(corr);
1008   corr->SetTitle("Trigger bias correction");
1009
1010   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
1011   InitPadCOLZ();
1012   corr->DrawCopy("COLZ");
1013
1014   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
1015   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
1016
1017   corr->GetYaxis()->SetRangeUser(0, 5);
1018
1019   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
1020   InitPadCOLZ();
1021   corr->DrawCopy("COLZ");
1022
1023   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
1024   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
1025 }
1026
1027 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1028 {
1029   gSystem->Load("libPWG0base");
1030
1031   TFile* file = TFile::Open(fileName);
1032   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1033   dNdEtaCorrection->LoadHistograms();
1034
1035   TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
1036   TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1037
1038   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
1039   canvas->Divide(2, 1);
1040
1041   canvas->cd(1);
1042   InitPad();
1043
1044   Prepare1DPlot(hist);
1045   hist->SetTitle("");
1046   hist->GetYaxis()->SetTitle("correction factor");
1047   hist->GetYaxis()->SetRangeUser(1, 1.5);
1048   hist->GetYaxis()->SetTitleOffset(1.6);
1049   hist->Draw();
1050
1051   canvas->cd(2);
1052   InitPad();
1053
1054   Prepare1DPlot(hist2);
1055   hist2->SetTitle("");
1056   hist2->GetYaxis()->SetTitle("correction factor");
1057   hist2->GetXaxis()->SetRangeUser(0, 5);
1058   hist2->GetYaxis()->SetTitleOffset(1.6);
1059   hist2->GetXaxis()->SetTitle("multiplicity");
1060   hist2->Draw();
1061
1062   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1063   pave->SetFillColor(0);
1064   pave->AddText("|z| < 10 cm");
1065   pave->Draw();
1066
1067   canvas->SaveAs("TriggerBias1D.eps");
1068 }
1069
1070 void VtxRecon()
1071 {
1072   TFile* file = TFile::Open("correction_map.root");
1073
1074   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
1075
1076   Prepare2DPlot(corr);
1077   corr->SetTitle("Vertex reconstruction correction");
1078
1079   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
1080   InitPadCOLZ();
1081   corr->DrawCopy("COLZ");
1082
1083   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
1084   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
1085
1086   corr->GetYaxis()->SetRangeUser(0, 5);
1087
1088   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
1089   InitPadCOLZ();
1090   corr->DrawCopy("COLZ");
1091
1092   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
1093   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
1094 }
1095
1096 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1097 {
1098   gSystem->Load("libPWG0base");
1099
1100   TFile* file = TFile::Open(fileName);
1101   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1102   dNdEtaCorrection->LoadHistograms();
1103
1104   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
1105   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1106
1107   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
1108   canvas->Divide(2, 1);
1109
1110   canvas->cd(1);
1111   InitPad();
1112
1113   Prepare1DPlot(hist);
1114   hist->SetTitle("");
1115   hist->GetYaxis()->SetTitle("correction factor");
1116   hist->GetYaxis()->SetRangeUser(1, 1.8);
1117   hist->GetYaxis()->SetTitleOffset(1.6);
1118   hist->DrawCopy();
1119
1120   canvas->cd(2);
1121   InitPad();
1122
1123   Prepare1DPlot(hist2);
1124   hist2->SetTitle("");
1125   hist2->GetYaxis()->SetTitle("correction factor");
1126   hist2->GetXaxis()->SetRangeUser(0, 20);
1127   hist2->GetYaxis()->SetTitleOffset(1.6);
1128   hist2->GetXaxis()->SetTitle("multiplicity");
1129   hist2->Draw();
1130
1131   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1132   pave->SetFillColor(0);
1133   pave->AddText("|z| < 10 cm");
1134   pave->Draw();
1135
1136   canvas->SaveAs("VtxRecon1D.eps");
1137
1138   Correction1DCreatePlots(fileName, folderName, 9.9, 2);
1139
1140   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1141   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1142
1143   Prepare1DPlot(corrX);
1144   Prepare1DPlot(corrZ);
1145
1146   corrX->GetYaxis()->SetTitleOffset(1.5);
1147   corrZ->GetYaxis()->SetTitleOffset(1.5);
1148
1149   corrX->SetTitle("a) z projection");
1150   corrZ->SetTitle("b) p_{T} projection");
1151
1152   corrX->GetYaxis()->SetTitle("Correction factor");
1153   corrZ->GetYaxis()->SetTitle("Correction factor");
1154
1155   corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
1156
1157   TString canvasName;
1158   canvasName.Form("VtxRecon1D_Track");
1159   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
1160   canvas->Divide(2, 1);
1161
1162   canvas->cd(1);
1163   InitPad();
1164   corrX->DrawCopy();
1165
1166   canvas->cd(2);
1167   InitPad();
1168   gPad->SetLogx();
1169   corrZ->Draw();
1170
1171   canvas->SaveAs("VtxRecon1D_Track.eps");
1172   canvas->SaveAs("VtxRecon1D_Track.gif");
1173 }
1174
1175 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
1176 {
1177   gSystem->Load("libPWG0base");
1178
1179   TFile::Open(fileName);
1180   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1181   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
1182
1183   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1184   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1185
1186   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1187   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1188
1189   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1190   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1191   gene->GetXaxis()->SetRangeUser(-10, 10);
1192   meas->GetXaxis()->SetRangeUser(-10, 10);
1193
1194   Float_t eff1 = gene->Integral() / meas->Integral();
1195   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1196
1197   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
1198
1199   gene->GetZaxis()->SetRangeUser(0.3, 10);
1200   meas->GetZaxis()->SetRangeUser(0.3, 10);
1201
1202   Float_t eff2 = gene->Integral() / meas->Integral();
1203   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1204
1205   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
1206
1207   gene->GetZaxis()->SetRangeUser(0.3, 1);
1208   meas->GetZaxis()->SetRangeUser(0.3, 1);
1209
1210   Float_t eff3 = gene->Integral() / meas->Integral();
1211   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1212
1213   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
1214 }
1215
1216 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)
1217 {
1218   if (correctionType2 == -1)
1219     correctionType2 = correctionType;
1220
1221   TFile::Open(fileName);
1222   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1223   dNdEtaCorrection->LoadHistograms();
1224
1225   TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
1226   TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType2)->GetTrackCorrection()->GetMeasuredHistogram();
1227
1228   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1229   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1230   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1231   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1232   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kFALSE);
1233   gene->GetYaxis()->SetRange(0, 0);
1234   meas->GetYaxis()->SetRange(0, 0);
1235
1236   gene->GetXaxis()->SetRangeUser(-9.9, 9.9);
1237   meas->GetXaxis()->SetRangeUser(-9.9, 9.9);
1238   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kFALSE);
1239   gene->GetZaxis()->SetRange(0, 0);
1240   meas->GetZaxis()->SetRange(0, 0);
1241
1242   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1243   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1244   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kFALSE);
1245 }
1246
1247 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)
1248 {
1249   gSystem->Load("libPWG0base");
1250
1251   Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType, correctionType2);
1252
1253   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1254   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1255   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1256
1257   Prepare1DPlot(corrX);
1258   Prepare1DPlot(corrY);
1259   Prepare1DPlot(corrZ);
1260
1261   /*
1262   corrX->SetTitle("a) z projection");
1263   corrY->SetTitle("b) #eta projection");
1264   corrZ->SetTitle("c) p_{T} projection");
1265   */
1266
1267   corrX->SetTitle("");
1268   corrY->SetTitle("");
1269   corrZ->SetTitle("");
1270
1271   corrX->SetTitleSize(0.06, "xyz");
1272   corrX->SetLabelSize(0.06, "xyz");
1273   corrY->SetTitleSize(0.06, "xyz");
1274   corrY->SetLabelSize(0.06, "xyz");
1275   corrZ->SetTitleSize(0.06, "xyz");
1276   corrZ->SetLabelSize(0.06, "xyz");
1277
1278   corrX->GetYaxis()->SetTitle("Correction factor");
1279   corrY->GetYaxis()->SetTitle("Correction factor");
1280   corrZ->GetYaxis()->SetTitle("Correction factor");
1281   //corrX->GetYaxis()->SetTitleOffset(1.7);
1282   //corrY->GetYaxis()->SetTitleOffset(1.7);
1283   //corrZ->GetYaxis()->SetTitleOffset(1.7);
1284   corrX->GetYaxis()->SetRangeUser(0.8, 1.5);
1285   corrY->GetYaxis()->SetRangeUser(0.8, 1.5);
1286   corrZ->GetYaxis()->SetRangeUser(0.8, 1.5);
1287
1288   corrZ->GetXaxis()->SetRangeUser(0.11, upperPtLimit);
1289
1290   TString canvasName;
1291   canvasName.Form(Form("Correction1D_%d_%s_%f", correctionType, fileName, upperPtLimit));
1292   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1293   canvas->Divide(3, 1);
1294
1295   TLatex* Tl = new TLatex;
1296   Tl->SetTextSize(0.06);
1297   Tl->SetBit(TLatex::kTextNDC);
1298
1299   canvas->cd(1);
1300   InitPad();
1301   gPad->SetTopMargin(0.05);
1302   gPad->SetBottomMargin(0.15);
1303   corrX->DrawCopy();
1304   Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1305   Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
1306
1307   canvas->cd(2);
1308   InitPad();
1309   gPad->SetTopMargin(0.05);
1310   gPad->SetBottomMargin(0.15);
1311   corrY->Draw();
1312   Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1313   Tl->DrawLatex(0.5, 0.8, "|vtx-z| < 10 cm");
1314
1315   canvas->cd(3);
1316   InitPad();
1317   gPad->SetTopMargin(0.05);
1318   gPad->SetBottomMargin(0.15);
1319   gPad->SetLogx();
1320   corrZ->Draw();
1321   corrZ->GetXaxis()->SetLabelOffset(0.005);
1322   corrZ->GetXaxis()->SetTitleOffset(1.2);
1323   Tl->DrawLatex(0.5, 0.88, "|vtx-z| < 10 cm");
1324   Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
1325
1326   return canvas;
1327 }
1328
1329 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1330 {
1331   gSystem->Load("libPWG0base");
1332
1333   Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
1334
1335   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1336   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1337   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1338
1339   Prepare1DPlot(corrX);
1340   Prepare1DPlot(corrY);
1341   Prepare1DPlot(corrZ);
1342
1343   corrX->SetTitle("a) z projection");
1344   corrY->SetTitle("a) #eta projection");
1345   corrZ->SetTitle("b) p_{T} projection");
1346
1347   corrY->GetYaxis()->SetTitle("correction factor");
1348   corrZ->GetYaxis()->SetTitle("correction factor");
1349
1350   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1351
1352   TString canvasName;
1353   canvasName.Form("Track2Particle1D_%s", folder);
1354   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1355   canvas->Divide(3, 1);
1356
1357   canvas->cd(1);
1358   InitPad();
1359   corrX->DrawCopy();
1360
1361   canvas->cd(2);
1362   InitPad();
1363   corrY->Draw();
1364
1365   canvas->cd(3);
1366   InitPad();
1367   corrZ->Draw();
1368
1369   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
1370   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
1371
1372   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
1373
1374   canvasName.Form("Track2Particle1D_%s_etapt", folder);
1375   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
1376   canvas->Divide(2, 1);
1377
1378   canvas->cd(1);
1379   InitPad();
1380   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
1381   corrY->GetYaxis()->SetRangeUser(1, 1.5);
1382   corrY->GetYaxis()->SetTitleOffset(1.5);
1383   corrY->DrawCopy();
1384   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
1385   pave->AddText("|z| < 10 cm");
1386   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
1387   pave->Draw();
1388
1389   canvas->cd(2);
1390   InitPad();
1391   gPad->SetLogx();
1392   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
1393   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
1394   corrZ->GetYaxis()->SetTitleOffset(1.5);
1395   corrZ->DrawCopy();
1396   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
1397   pave->AddText("|z| < 10 cm");
1398   pave->AddText("|#eta| < 0.8");
1399   pave->Draw();
1400
1401   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
1402   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
1403 }
1404
1405 /*
1406 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1407 {
1408   gSystem->Load("libPWG0base");
1409
1410   // particle type
1411   for (Int_t particle=0; particle<4; ++particle)
1412   {
1413     TString dirName;
1414     dirName.Form("correction_%d", particle);
1415     Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
1416
1417     TString tmpx, tmpy, tmpz;
1418     tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1419     tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1420     tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
1421
1422     TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1423     TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1424     TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
1425
1426     Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
1427
1428     TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1429     TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1430     TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
1431
1432     posX->Divide(negX);
1433     posY->Divide(negY);
1434     posZ->Divide(negZ);
1435
1436     Prepare1DPlot(posX);
1437     Prepare1DPlot(posY);
1438     Prepare1DPlot(posZ);
1439
1440     Float_t min = 0.8;
1441     Float_t max = 1.2;
1442
1443     posX->SetMinimum(min);
1444     posX->SetMaximum(max);
1445     posY->SetMinimum(min);
1446     posY->SetMaximum(max);
1447     posZ->SetMinimum(min);
1448     posZ->SetMaximum(max);
1449
1450     posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1451
1452     posX->GetYaxis()->SetTitleOffset(1.7);
1453     posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1454     posY->GetYaxis()->SetTitleOffset(1.7);
1455     posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1456     posZ->GetYaxis()->SetTitleOffset(1.7);
1457     posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
1458
1459     posZ->GetXaxis()->SetRangeUser(0, 1);
1460
1461     TString canvasName;
1462     canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
1463
1464     TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1465     canvas->Divide(3, 1);
1466
1467     canvas->cd(1);
1468     InitPad();
1469     posX->DrawCopy();
1470
1471     canvas->cd(2);
1472     InitPad();
1473     posY->DrawCopy();
1474
1475     canvas->cd(3);
1476     InitPad();
1477     posZ->DrawCopy();
1478
1479     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1480     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1481   }
1482 }
1483 */
1484
1485 void CompareTrack2Particle1D(const char* file1, const char* file2, Float_t upperPtLimit = 9.9)
1486 {
1487   loadlibs();
1488
1489   const char* folderName = "dndeta_correction";
1490
1491   c = new TCanvas("CompareTrack2Particle1D", "CompareTrack2Particle1D", 1200, 400);
1492   c->Divide(3, 1);
1493
1494   for (Int_t fileId = 0; fileId < 2; fileId++)
1495   {
1496     const char* file = ((fileId == 0) ? file1 : file2);
1497     Correction1DCreatePlots(file, folderName, upperPtLimit, 1);
1498
1499     TH1* corr[3];
1500     corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1501     corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y_div_measured_y"));
1502     corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1503     /*corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x"))->Clone(Form("hist_x_%d", fileId));
1504     corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y"))->Clone(Form("hist_y_%d", fileId));
1505     corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z"))->Clone(Form("hist_z_%d", fileId));*/
1506
1507     for (Int_t i=0; i<3; i++)
1508     {
1509       c->cd(i+1);
1510       InitPad();
1511       corr[i]->GetYaxis()->SetRangeUser(0.8, 2);
1512       corr[i]->SetLineColor(fileId+1);
1513       corr[i]->DrawCopy((fileId == 0) ? "" : "SAME");
1514     }
1515   }
1516
1517   return;
1518
1519   c->SaveAs(Form("%s.gif", canvas->GetName()));
1520   c->SaveAs(Form("%s.eps", canvas->GetName()));
1521 }
1522
1523 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
1524 {
1525   TFile::Open(fileName);
1526   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1527   dNdEtaCorrection->LoadHistograms();
1528
1529   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1530   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1531
1532   gene->GetZaxis()->SetRangeUser(0.2, 10);
1533   meas->GetZaxis()->SetRangeUser(0.2, 10);
1534   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
1535   gene->GetZaxis()->SetRange(0, 0);
1536   meas->GetZaxis()->SetRange(0, 0);
1537
1538   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1539   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1540   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
1541   gene->GetYaxis()->SetRange(0, 0);
1542   meas->GetYaxis()->SetRange(0, 0);
1543
1544   gene->GetXaxis()->SetRangeUser(-10, 10);
1545   meas->GetXaxis()->SetRangeUser(-10, 10);
1546   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
1547   gene->GetXaxis()->SetRange(0, 0);
1548   meas->GetXaxis()->SetRange(0, 0);
1549 }
1550
1551 TCanvas* Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1552 {
1553   gSystem->Load("libPWG0base");
1554
1555   Track2Particle2DCreatePlots(fileName);
1556
1557   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1558   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1559   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
1560
1561   Prepare2DPlot(corrYX);
1562   Prepare2DPlot(corrZX);
1563   Prepare2DPlot(corrZY);
1564
1565   const char* title = "";
1566   corrYX->SetTitle(title);
1567   corrZX->SetTitle(title);
1568   corrZY->SetTitle(title);
1569
1570   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1571   canvas->Divide(3, 1);
1572
1573   canvas->cd(1);
1574   InitPadCOLZ();
1575   corrYX->Draw("COLZ");
1576
1577   canvas->cd(2);
1578   InitPadCOLZ();
1579   corrZX->Draw("COLZ");
1580
1581   canvas->cd(3);
1582   InitPadCOLZ();
1583   corrZY->Draw("COLZ");
1584
1585   canvas->SaveAs(Form("corr_track2particle_%d.gif", gMax));
1586   canvas->SaveAs(Form("corr_track2particle_%d.eps", gMax));
1587
1588   return canvas;
1589 }
1590
1591 void CompareTrack2Particle2D()
1592 {
1593   gSystem->Load("libPWG0base");
1594
1595   Track2Particle2DCreatePlots("correction_maponly-positive.root");
1596
1597   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("pos_yx"));
1598   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("pos_zx"));
1599   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("pos_zy"));
1600
1601   Track2Particle2DCreatePlots("correction_maponly-negative.root");
1602
1603   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("neg_yx"));
1604   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("neg_zx"));
1605   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("neg_zy"));
1606
1607   posYX->Divide(negYX);
1608   posZX->Divide(negZX);
1609   posZY->Divide(negZY);
1610
1611   Prepare2DPlot(posYX);
1612   Prepare2DPlot(posZX);
1613   Prepare2DPlot(posZY);
1614
1615   Float_t min = 0.8;
1616   Float_t max = 1.2;
1617
1618   posYX->SetMinimum(min);
1619   posYX->SetMaximum(max);
1620   posZX->SetMinimum(min);
1621   posZX->SetMaximum(max);
1622   posZY->SetMinimum(min);
1623   posZY->SetMaximum(max);
1624
1625   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1626   canvas->Divide(3, 1);
1627
1628   canvas->cd(1);
1629   InitPadCOLZ();
1630   posYX->Draw("COLZ");
1631
1632   canvas->cd(2);
1633   InitPadCOLZ();
1634   posZX->Draw("COLZ");
1635
1636   canvas->cd(3);
1637   InitPadCOLZ();
1638   posZY->Draw("COLZ");
1639
1640   canvas->SaveAs("CompareTrack2Particle2D.gif");
1641   canvas->SaveAs("CompareTrack2Particle2D.eps");
1642 }
1643
1644 void Track2Particle3D()
1645 {
1646   // get left margin proper
1647
1648   TFile* file = TFile::Open("correction_map.root");
1649
1650   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1651
1652   corr->SetTitle("Correction Factor");
1653   SetRanges(corr->GetZaxis());
1654
1655   Prepare3DPlot(corr);
1656
1657   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1658   canvas->SetTheta(29.428);
1659   canvas->SetPhi(16.5726);
1660
1661   corr->Draw();
1662
1663   canvas->SaveAs("Track2Particle3D.gif");
1664   canvas->SaveAs("Track2Particle3D.eps");
1665 }
1666
1667 void Track2Particle3DAll()
1668 {
1669   TFile* file = TFile::Open("correction_map.root");
1670
1671   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1672   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1673   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1674
1675   gene->SetTitle("Generated Particles");
1676   meas->SetTitle("Measured Tracks");
1677   corr->SetTitle("Correction Factor");
1678
1679   Prepare3DPlot(gene);
1680   Prepare3DPlot(meas);
1681   Prepare3DPlot(corr);
1682
1683   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1684   canvas->Divide(3, 1);
1685
1686   canvas->cd(1);
1687   InitPad();
1688   gene->Draw();
1689
1690   canvas->cd(2);
1691   meas->Draw();
1692
1693   canvas->cd(3);
1694   corr->Draw();
1695
1696   canvas->SaveAs("Track2Particle3DAll.gif");
1697   canvas->SaveAs("Track2Particle3DAll.eps");
1698 }
1699
1700 void MultiplicityMC(Int_t xRangeMax = 50)
1701 {
1702   TFile* file = TFile::Open("multiplicityMC.root");
1703
1704   if (!file)
1705   {
1706     printf("multiplicityMC.root could not be opened.\n");
1707     return;
1708   }
1709
1710   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1711   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1712   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1713
1714   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1715   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1716   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1717   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1718   {
1719     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1720     proj->Fit("gaus", "0");
1721     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1722     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1723
1724     continue;
1725
1726     // draw for debugging
1727     new TCanvas;
1728     proj->DrawCopy();
1729     proj->GetFunction("gaus")->DrawCopy("SAME");
1730   }
1731
1732   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1733
1734   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1735   {
1736     Float_t mean = correction->GetBinContent(i);
1737     Float_t width = correctionWidth->GetBinContent(i);
1738
1739     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1740     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1741     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1742
1743     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1744     {
1745       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1746     }
1747   }
1748
1749   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1750   fMultiplicityESDCorrectedRebinned->Rebin(10);
1751   fMultiplicityESDCorrectedRebinned->Scale(0.1);
1752
1753   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1754   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1755   ratio->Divide(fMultiplicityMC);
1756
1757   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1758   ratio2->Divide(fMultiplicityMC);
1759
1760   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1761   canvas->Divide(3, 2);
1762
1763   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1764   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1765   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1766   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1767   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1768   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1769   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1770
1771   canvas->cd(1); //InitPad();
1772   fMultiplicityESD->Draw();
1773   fMultiplicityMC->SetLineColor(2);
1774   fMultiplicityMC->Draw("SAME");
1775
1776   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1777   legend->AddEntry(fMultiplicityESD, "ESD");
1778   legend->AddEntry(fMultiplicityMC, "MC");
1779   legend->Draw();
1780
1781   canvas->cd(2);
1782   fCorrelation->Draw("COLZ");
1783
1784   canvas->cd(3);
1785   correction->Draw();
1786   //correction->Fit("pol1");
1787   correctionWidth->SetLineColor(2);
1788   correctionWidth->Draw("SAME");
1789
1790   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1791   legend->AddEntry(correction, "#bar{x}");
1792   legend->AddEntry(correctionWidth, "#sigma");
1793   legend->Draw();
1794
1795   canvas->cd(4);
1796   ratio->Draw();
1797
1798   ratio2->SetLineColor(2);
1799   ratio2->Draw("SAME");
1800
1801   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1802   legend->AddEntry(ratio, "uncorrected");
1803   legend->AddEntry(ratio2, "corrected");
1804   legend->Draw();
1805
1806   canvas->cd(5);
1807   fMultiplicityESDCorrected->SetLineColor(kBlue);
1808   fMultiplicityESDCorrected->Draw();
1809   fMultiplicityMC->Draw("SAME");
1810   fMultiplicityESD->Draw("SAME");
1811
1812   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1813   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1814   legend->AddEntry(fMultiplicityMC, "MC");
1815   legend->AddEntry(fMultiplicityESD, "ESD");
1816   legend->Draw();
1817
1818   canvas->cd(6);
1819   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
1820   fMultiplicityESDCorrectedRebinned->Draw();
1821   fMultiplicityMC->Draw("SAME");
1822
1823   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1824   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1825   legend->AddEntry(fMultiplicityMC, "MC");
1826   legend->Draw();
1827
1828   canvas->SaveAs("MultiplicityMC.gif");
1829 }
1830
1831 void MultiplicityESD()
1832 {
1833   TFile* file = TFile::Open("multiplicityESD.root");
1834
1835   if (!file)
1836   {
1837     printf("multiplicityESD.root could not be opened.\n");
1838     return;
1839   }
1840
1841   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1842
1843   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1844
1845   fMultiplicityESD->Draw();
1846 }
1847
1848 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")
1849 {
1850   loadlibs();
1851
1852   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1853   TFile::Open(correctionMapFile);
1854   dNdEtaCorrection->LoadHistograms();
1855
1856   TFile* file = TFile::Open(dataInput);
1857
1858   if (!file)
1859   {
1860     cout << "Error. File not found" << endl;
1861     return;
1862   }
1863
1864   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1865   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1866
1867   gROOT->cd();
1868
1869   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
1870   hist1->SetTitle("mc");
1871   Printf("mc contains %f entries", hist1->Integral());
1872   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()));
1873
1874   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
1875   hist2->SetTitle("esd");
1876   Printf("esd contains %f entries", hist2->Integral());
1877   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()));
1878
1879   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1880   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
1881
1882   hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1883   hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1884   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
1885
1886   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
1887   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
1888   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
1889   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
1890   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
1891
1892   TH2* hist3 = (TH2*) dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram()->Clone("mc2");
1893   hist3->SetTitle("mc2");
1894   Printf("mc event contains %f entries", hist3->Integral());
1895   Printf("mc event contains %f entries in |vtx-z| < 10", hist3->Integral(hist3->GetXaxis()->FindBin(-9.9), hist3->GetXaxis()->FindBin(9.9), 1, hist3->GetNbinsY()));
1896
1897   TH2* hist4 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("esd2");
1898   hist4->SetTitle("esd2");
1899   Printf("esd event contains %f entries", hist4->Integral());
1900   Printf("esd event contains %f entries in |vtx-z| < 10", hist4->Integral(hist4->GetXaxis()->FindBin(-9.9), hist4->GetXaxis()->FindBin(9.9), 1, hist4->GetNbinsY()));
1901
1902   ratio = (TH2*) hist3->Clone("ratio");
1903   ratio->Divide(hist4);
1904
1905   new TCanvas; ratio->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 void DrawTrackletOrigin_Compare(const char* file1, const char* file2)
2226 {
2227   DrawTrackletOrigin(file1);
2228   good1 = (TH1*) gROOT->FindObject("good")->Clone("good1");
2229   bad1 = (TH1*) gROOT->FindObject("bad")->Clone("bad1");
2230
2231   DrawTrackletOrigin(file2);
2232   good2 = (TH1*) gROOT->FindObject("good")->Clone("good2");
2233   bad2 = (TH1*) gROOT->FindObject("bad")->Clone("bad2");
2234
2235   c = new TCanvas("c4", "c4", 600, 600);
2236   c->SetRightMargin(0.05);
2237   c->SetTopMargin(0.05);
2238   c->SetGridx();
2239   c->SetGridy();
2240   gPad->SetLogy();
2241
2242   good1->Draw();
2243   bad1->SetLineColor(1);
2244   bad1->SetMarkerColor(1);
2245   bad1->Draw("SAME");
2246
2247   Float_t factor = (good1->Integral() + bad1->Integral()) / (good2->Integral() + bad2->Integral());
2248
2249   good2->Scale(factor);
2250   bad2->Scale(factor);
2251
2252   good2->SetLineColor(2);
2253   bad2->SetMarkerColor(2);
2254
2255   good2->Draw("SAME");
2256   bad2->Draw("SAME");
2257
2258   good1->GetYaxis()->SetRangeUser(1, TMath::Max(good1->GetMaximum(), good2->GetMaximum()) * 1.1);
2259 }
2260
2261 void Tracklets_Asymmetry()
2262 {
2263   TFile::Open("correction_map.root");
2264
2265   Int_t maxHists = 7;
2266   TH1* hist[8];
2267
2268   Int_t colors[]  = {1,2,3,4,6,7,8,102};
2269   const char* titles[] = { "PP", "SS", "PP'", "PS'", "PS", "SP'", "SS'", "" };
2270
2271   TLegend* legend = new TLegend(0.75, 0.6, 0.93, 0.93);
2272
2273   for (Int_t i=0; i<maxHists; i++)
2274   {
2275     hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
2276     hist[i]->Rebin(10);
2277
2278     for (Int_t j=hist[i]->GetNbinsX()/2; j<=hist[i]->GetNbinsX(); j++)
2279       if (hist[i]->GetBinContent(j) > 0)
2280         hist[i]->SetBinContent(j, (hist[i]->GetBinContent(j) -  hist[i]->GetBinContent(hist[i]->GetXaxis()->FindBin(-hist[i]->GetXaxis()->GetBinCenter(j)))) / hist[i]->GetBinContent(j));
2281
2282     hist[i]->SetStats(kFALSE);
2283     hist[i]->SetLineColor(colors[i]);
2284     hist[i]->GetXaxis()->SetRangeUser(0.001, 0.09);
2285     //hist[i]->GetYaxis()->SetRangeUser(5, hist[i]->GetMaximum() * 2);
2286     hist[i]->GetYaxis()->SetTitleOffset(1.3);
2287     hist[i]->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2288     hist[i]->Draw(((i == 0) ? "" : "SAME"));
2289
2290     legend->AddEntry(hist[i], titles[i], "L");
2291   }
2292
2293   legend->Draw();
2294 }
2295
2296 TH2* GetCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2297 {
2298   // returns the correction factor with pt integrated out
2299
2300   loadlibs();
2301
2302   TFile::Open(fileName);
2303
2304   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
2305   if (!dNdEtaCorrection->LoadHistograms())
2306     return;
2307
2308   //  hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
2309
2310   gener = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
2311   measu = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
2312
2313   gener->GetZaxis()->SetRange(gener->GetZaxis()->FindBin(ptmin), gener->GetNbinsZ()+1);
2314   TH2D *gener_xy = gener->Project3D("yx");
2315
2316   measu->GetZaxis()->SetRange(measu->GetZaxis()->FindBin(ptmin), measu->GetNbinsZ()+1);
2317   TH2D *measu_xy = measu->Project3D("yx");
2318
2319   cout << measu->GetZaxis()->FindBin(ptmin) << " " << measu->GetNbinsZ()+1 << endl;
2320
2321   TCanvas *canp = new TCanvas("canp","canp",600,1000);
2322   canp->Divide(1,2,0.0001,0.0001);
2323   canp->cd(1);
2324   gener_xy->Draw("COLZ");
2325   canp->cd(2);
2326   measu_xy->Draw("COLZ");
2327
2328
2329   TCanvas *canpr = new TCanvas("canpr","canpr",700,500);
2330   canpr->cd();
2331   TH2D *proj = new TH2D(*gener_xy);
2332   proj->Divide(measu_xy);
2333
2334 //   proj = hist->Project3D("yx");
2335   proj->Draw("COLZ");
2336
2337   return proj;
2338 }
2339
2340 void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2341 {
2342   TH2* proj = GetCorrection(fileName, dirName, ptmin);
2343
2344   const Float_t limit = 5;
2345
2346   TString array = "{";
2347   TString arrayEnd = "}";
2348
2349   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2350   {
2351     Int_t begin = -1;
2352     Int_t end = -1;
2353     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2354     {
2355       if (begin == -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2356         begin = x;
2357       if (begin != -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2358         end = x;
2359     }
2360     Printf("Limits for y = %d are %d to %d", y, begin, end);
2361
2362     if (y > 1)
2363       array += ", ";
2364     array += Form("%d", begin);
2365
2366     if (y > 1)
2367       arrayEnd.Prepend(", ");
2368     arrayEnd.Prepend(Form("%d", (end == -1) ? -1 : proj->GetNbinsX() + 1 - end));
2369   }
2370   array += "}";
2371   arrayEnd.Prepend("{");
2372
2373   Printf("Begin array:");
2374   Printf("%s", array.Data());
2375
2376   Printf("End array (mirrored) (should be the same):");
2377   Printf("%s", arrayEnd.Data());
2378 }
2379
2380 void AverageMultiplicity(const char* fileName = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
2381 {
2382   loadlibs();
2383
2384   TFile::Open(fileName);
2385
2386   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
2387   dNdEtaCorrection->LoadHistograms();
2388   TH2* events = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram();
2389   TH3* tracks = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram();
2390
2391   Float_t nEvents = events->Integral(events->GetXaxis()->FindBin(-1), events->GetXaxis()->FindBin(1), 0, events->GetNbinsY()+1);
2392   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);
2393
2394   Printf("%f %f --> %f", nEvents, nTracks, nTracks / nEvents);
2395 }
2396
2397 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")
2398 {
2399   loadlibs();
2400
2401   TFile::Open(rawFile);
2402   dNdEtaAnalysis* raw = new dNdEtaAnalysis("dndeta", "dndeta");
2403   raw->LoadHistograms("fdNdEtaAnalysisESD");
2404   raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->GetXaxis()->SetRangeUser(-vertexRange, vertexRange);
2405   tracks = raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Project3D("y");
2406   events = raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("events", 0, raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->GetNbinsY() + 1);
2407   Float_t nEvents = events->Integral(events->FindBin(-vertexRange), events->FindBin(vertexRange));
2408   tracks->Scale(1.0 / nEvents / tracks->GetBinWidth(1));
2409
2410   TFile::Open(mcFile);
2411   dNdEtaAnalysis* mc = new dNdEtaAnalysis("dndeta", "dndeta");
2412   mc->LoadHistograms("dndetaTrVtx");
2413   mcH = mc->GetdNdEtaPtCutOffCorrectedHistogram(0);
2414
2415   new TCanvas;
2416   mcH->SetLineColor(2);
2417   mcH->DrawCopy();
2418   tracks->DrawCopy("SAME");
2419
2420   new TCanvas;
2421   mcH->GetYaxis()->SetRangeUser(0, 5);
2422   mcH->Divide(tracks);
2423   mcH->DrawCopy();
2424   mcH->Fit("pol0", "", "", -etaRange, etaRange);
2425 }
2426
2427 void TrackCuts_Comparison_MC(char* histName, Int_t plotWhich = 0, const char* fileName = "correction_map.root", Bool_t mirror = kFALSE)
2428 {
2429   // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
2430   //    --> manually disable it in the run.C
2431   //
2432   // plotWhich: 0 = only before
2433   //            1 = both
2434   //            2 = only after
2435   //
2436   // mirror: kTRUE --> project negative values on the positive side
2437
2438
2439   file = TFile::Open(fileName);
2440
2441   Int_t count = 0;
2442   Int_t colors[] = { 1, 2, 3, 4, 5, 6 };
2443
2444   TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
2445   TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
2446   TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
2447
2448   TCanvas* c1 = new TCanvas(histName, histName, 800, 1200);
2449   c1->Divide(1, 2);
2450   //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
2451   //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
2452
2453   const char* folders2[] = { "before_cuts", "after_cuts" };
2454   Bool_t first = kTRUE;
2455   for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
2456   {
2457     const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
2458     const char* names[] =    { "all", "primaries", "secondaries" };
2459     TH1* base = 0;
2460     TH1* prim = 0;
2461     TH1* sec = 0;
2462     for (Int_t i = 0; i < 3; i++)
2463     {
2464       TString folder;
2465       folder.Form("%s/%s/%s", folders1[i], folders2[j], histName);
2466       TH1* hist = (TH1*) file->Get(folder);
2467
2468       if (mirror)
2469       {
2470         for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
2471         {
2472           Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
2473           if (bin != newBin)
2474           {
2475             hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
2476             hist->SetBinContent(bin, 0);
2477           }
2478         }
2479       }
2480
2481       legend->AddEntry(hist, Form("%s %s", names[i], folders2[j]));
2482
2483       c1->cd(1);
2484       hist->SetLineColor(colors[count]);
2485       hist->DrawCopy((count == 0) ? "" : "SAME");
2486
2487       switch (i)
2488       {
2489         case 0: base = hist; break;
2490         case 1: prim = hist; break;
2491         case 2: sec = hist; break;
2492       }
2493
2494       count++;
2495     }
2496
2497     TH1* eff    = (TH1*) prim->Clone("eff"); eff->Reset();
2498     TH1* purity = (TH1*) prim->Clone("purity"); purity->Reset();
2499
2500     for (Int_t bin = 1; bin <= prim->GetNbinsX(); bin++)
2501     {
2502       eff->SetBinContent(bin, prim->Integral(1, bin) / prim->Integral(1, prim->GetNbinsX() + 1));
2503       if (prim->Integral(1, bin) + sec->Integral(1, bin) > 0)
2504         purity->SetBinContent(bin, sec->Integral(1, bin) / (prim->Integral(1, bin) + sec->Integral(1, bin)));
2505     }
2506
2507     eff->GetYaxis()->SetRangeUser(0, 1);
2508     eff->SetLineColor(colors[0+j*2]);
2509     eff->SetStats(kFALSE);
2510     purity->SetLineColor(colors[1+j*2]);
2511
2512     legend3->AddEntry(eff, Form("%s: efficiency", folders2[j]));
2513     legend3->AddEntry(purity, Form("%s: contamination", folders2[j]));
2514
2515     c1->cd(2);
2516     eff->DrawCopy((first) ? "" : "SAME");
2517     first = kFALSE;
2518     purity->DrawCopy("SAME");
2519   }
2520
2521   c1->cd(1)->SetLogy();
2522   c1->cd(1)->SetGridx();
2523   c1->cd(1)->SetGridy();
2524   legend->Draw();
2525
2526   //c2->cd();
2527  // c2->SetGridx();
2528  // c2->SetGridy();
2529   //legend2->Draw();
2530
2531   c1->cd(2)->SetGridx();
2532   c1->cd(2)->SetGridy();
2533   legend3->Draw();
2534
2535   //c1->SaveAs(Form("%s.png", histName));
2536 }
2537
2538 void TrackCuts_Comparison_Data(char* histName, Int_t plotWhich, const char* fileName1, const char* fileName2, Bool_t mirror = kFALSE, const char* label1 = "file1", const char* label2 = "file2")
2539 {
2540   // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
2541   //    --> manually disable it in the run.C
2542   //
2543   // plotWhich: 0 = only before
2544   //            1 = both
2545   //            2 = only after
2546   //
2547   // mirror: kTRUE --> project negative values on the positive side
2548
2549
2550   Int_t count = 0;
2551   Int_t colors[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
2552
2553   TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
2554   legend->SetFillColor(0);
2555   TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
2556   TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
2557
2558   TCanvas* c1 = new TCanvas(histName, histName, 600, 600);
2559   //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
2560   //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
2561
2562   const char* folders2[] = { "before_cuts", "after_cuts" };
2563   Bool_t first = kTRUE;
2564   for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
2565   {
2566     const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
2567     const char* names[] =    { "all", "primaries", "secondaries" };
2568
2569     Float_t normalize[3];
2570
2571     for (Int_t i = 0; i < 2; i++)
2572     {
2573       file = TFile::Open((i == 0) ? fileName1 : fileName2);
2574
2575       for (Int_t k = 1; k < 3; k++)
2576       {
2577         TString folder;
2578         folder.Form("%s/%s/%s", folders1[k], folders2[j], histName);
2579         Printf("%s", folder.Data());
2580         TH1* hist = (TH1*) file->Get(folder);
2581
2582         if (mirror)
2583         {
2584           for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
2585           {
2586             Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
2587             if (bin != newBin)
2588             {
2589               hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
2590               hist->SetBinContent(bin, 0);
2591             }
2592           }
2593         }
2594
2595         if (i == 0)
2596         {
2597           normalize[k] = hist->Integral();
2598         }
2599         else
2600           hist->Scale(normalize[k] / hist->Integral());
2601
2602         legend->AddEntry(hist, Form("%s %s %s", (i == 0) ? label1 : label2, (k == 1) ? "primaries" : "secondaries", folders2[j]));
2603
2604         c1->cd();
2605         hist->SetStats(0);
2606         hist->SetLineColor(colors[count]);
2607         hist->DrawCopy((count == 0) ? "" : "SAME");
2608
2609         count++;
2610       }
2611     }
2612
2613   }
2614
2615   //c1->SetLogy();
2616   c1->SetGridx();
2617   c1->SetGridy();
2618   legend->Draw();
2619 }
2620
2621 void TrackCuts_DCA()
2622 {
2623   file = TFile::Open("correction_map.root");
2624   hist = (TH2*) file->Get("esd_track_cuts/before_cuts/dXYvsDZ");
2625
2626   TCanvas* c1 = new TCanvas("c1", "c1", 600, 600);
2627   c1->SetLogz();
2628   c1->SetRightMargin(0.12);
2629   c1->SetBottomMargin(0.12);
2630
2631   hist->SetStats(kFALSE);
2632   hist->Draw("COLZ");
2633
2634   ellipse = new TEllipse(0, 0, 4);
2635   ellipse->SetLineWidth(2);
2636   ellipse->SetLineStyle(2);
2637   ellipse->SetFillStyle(0);
2638   ellipse->Draw();
2639
2640   c1->SaveAs("trackcuts_dca_2d.eps");
2641 }
2642
2643 void FindNSigma(TH2* hist, Int_t nSigma = 1)
2644 {
2645   TH1* proj = hist->ProjectionY();
2646   proj->Reset();
2647
2648   for (Int_t bin=1; bin<=proj->GetNbinsX(); bin++)
2649   {
2650     if (hist->Integral(1, hist->GetNbinsX(), bin, bin) == 0)
2651       continue;
2652
2653     Int_t limit = -1;
2654     for (limit = 1; limit<=hist->GetNbinsX(); limit++)
2655   }
2656 }
2657
2658 void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2659 {
2660   TH2* proj = GetCorrection(fileName, dirName, ptmin);
2661
2662   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2663     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2664       if (proj->GetBinContent(x, y) > 5 || proj->GetBinContent(x, y) == 0)
2665       {
2666         proj->SetBinContent(x, y, 0);
2667       }
2668       else
2669         proj->SetBinContent(x, y, 1);
2670
2671
2672   input->Multiply(proj);
2673 }
2674
2675 void MakeGaussianProfile(const char* histName = "fVertexCorrelation", Bool_t subtractMean = kFALSE)
2676 {
2677     TFile::Open("correction_map.root");
2678
2679     TH2* hist2d = (TH2*) gFile->Get(histName);
2680     hist2d->Sumw2();
2681
2682     TH1* result = hist2d->ProjectionX("result");
2683     result->GetYaxis()->SetTitle(hist2d->GetYaxis()->GetTitle());
2684     result->Reset();
2685
2686     for (Int_t x=1; x<hist2d->GetNbinsX(); ++x)
2687     {
2688         hist = hist2d->ProjectionY(Form("temp_%d", x), x, x);
2689         if (hist->GetEntries() == 0)
2690             continue;
2691         if (hist->Fit("gaus") == 0)
2692         {
2693             func = hist->GetFunction("gaus");
2694             mean = func->GetParameter(1);
2695             error = func->GetParError(1);
2696
2697             if (subtractMean)
2698                 mean = hist2d->GetXaxis()->GetBinCenter(x) - mean;
2699
2700             result->SetBinContent(x, mean);
2701             result->SetBinError(x, error);
2702
2703             if (x % 10 == 0)
2704             {
2705                 new TCanvas;
2706                 ((TH1*) hist->Clone())->DrawCopy();
2707             }
2708         }
2709         //break;
2710     }
2711
2712     new TCanvas;
2713     result->GetYaxis()->SetRangeUser(-0.2, 0.2);
2714     result->Draw();
2715 }
2716
2717 TH2* GetAcceptance(void* corr2d_void)
2718 {
2719         corr2d = (AliCorrectionMatrix2D*) corr2d_void;
2720         corr_xy = (TH2*) corr2d->GetCorrectionHistogram()->Clone("acceptance");
2721
2722         // fold in acceptance
2723         for (Int_t x=1; x<=corr_xy->GetNbinsX(); ++x)
2724                 for (Int_t y=1; y<=corr_xy->GetNbinsY(); ++y)
2725                 {
2726                         if (corr_xy->GetBinContent(x, y) > 1.5)
2727                                 corr_xy->SetBinContent(x, y, 0);
2728
2729                         if (corr_xy->GetBinContent(x, y) > 0)
2730                                 corr_xy->SetBinContent(x, y, 1);
2731
2732                         corr_xy->SetBinError(x, y, 0);
2733                 }
2734
2735         return corr_xy;
2736 }
2737
2738 void ZeroOutsideAcceptance(TH2* acc, TH3* hist)
2739 {
2740   for (Int_t x=0; x<=acc->GetNbinsX()+1; ++x)
2741     for (Int_t y=0; y<=acc->GetNbinsY()+1; ++y)
2742     {
2743       if (acc->GetBinContent(x, y) > 1.5 || acc->GetBinContent(x, y) == 0)
2744       {
2745         for (Int_t z=0; z<=hist->GetNbinsZ()+1; ++z)
2746         {
2747           hist->SetBinContent(x, y, z, 0);
2748           hist->SetBinError(x, y, z, 0);
2749         }
2750       }
2751     }
2752 }
2753
2754 void DrawPhi()
2755 {
2756   loadlibs();
2757
2758   TFile::Open("correction_map.root");
2759   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
2760   if (!dNdEtaCorrection->LoadHistograms())
2761     return 0;
2762
2763   TFile::Open("analysis_esd.root");
2764   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
2765   fdNdEtaAnalysis->LoadHistograms();
2766
2767   // acc. map!
2768   //acc = GetAcceptance(dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000));
2769   acc = dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000)->GetCorrectionHistogram();
2770   //new TCanvas; acc->Draw("COLZ");
2771
2772   histG = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
2773   ZeroOutsideAcceptance(acc, histG);
2774   //new TCanvas; histG->Project3D("yx")->Draw("COLZ");
2775   //histG->GetYaxis()->SetRangeUser(-0.9, 0.9);
2776
2777   histM = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
2778   ZeroOutsideAcceptance(acc, histM);
2779   //histM->GetYaxis()->SetRangeUser(-0.9, 0.9);
2780
2781   TFile::Open("analysis_mc.root");
2782   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndetaTrVtxMC", "dndetaTrVtxMC");
2783   fdNdEtaAnalysis2->LoadHistograms("dndetaTrVtx");
2784
2785   histMC = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
2786   ZeroOutsideAcceptance(acc, histMC);
2787   //new TCanvas; histMC->Project3D("yx2")->Draw("COLZ");
2788
2789   //histG->GetZaxis()->SetRangeUser(1,2); histMC->GetZaxis()->SetRangeUser(1,2);
2790   new TCanvas; a = histG->Project3D("yx3"); a->Add(histMC->Project3D("yx4"), -1); a->Draw("COLZ");
2791
2792   //histMC->GetYaxis()->SetRangeUser(-0.9, 0.9);
2793
2794   c = new TCanvas;
2795
2796   histG->GetXaxis()->SetRangeUser(-9.9, 9.9);
2797   histG->Project3D("z")->Draw();
2798
2799   histM->GetXaxis()->SetRangeUser(-9.9, 9.9);
2800   proj = histM->Project3D("z2");
2801   proj->SetLineColor(2);
2802   proj->Draw("SAME");
2803
2804   histMC->GetXaxis()->SetRangeUser(-9.9, 9.9);
2805   projMC = histMC->Project3D("z3");
2806   projMC->SetLineColor(3);
2807   projMC->Draw("SAME");
2808 }
2809
2810 void PrintEventStats(Int_t corrID = 3)
2811 {
2812   loadlibs();
2813
2814   /*
2815   TFile::Open("analysis_mc.root");
2816   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD");
2817   fdNdEtaAnalysis->LoadHistograms();
2818   trackHist = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
2819   eventHist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetGeneratedHistogram();
2820   */
2821
2822   TFile::Open("correction_map.root");
2823   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
2824   if (!dNdEtaCorrection->LoadHistograms())
2825     return;
2826   trackHist = dNdEtaCorrection->GetCorrection(corrID)->GetTrackCorrection()->GetGeneratedHistogram();
2827   eventHist = dNdEtaCorrection->GetCorrection(corrID)->GetEventCorrection()->GetGeneratedHistogram();
2828
2829   trackHist->GetXaxis()->SetRange(0, trackHist->GetNbinsX()+1);
2830   trackHist->GetZaxis()->SetRange(0, trackHist->GetNbinsZ()+1);
2831   eta = trackHist->Project3D("y");
2832
2833   events = eventHist->Integral(0, eventHist->GetNbinsX()+1, 0, eventHist->GetNbinsY()+1);
2834
2835   eta->Scale(1.0 / events);
2836
2837   Float_t avgN = eta->Integral(0, eta->GetNbinsX()+1);
2838   Printf("<N> = %f", avgN);
2839
2840   eta->Scale(1.0 / eta->GetXaxis()->GetBinWidth(1));
2841
2842   Printf("dndeta | eta = 0 is %f", (eta->GetBinContent(eta->FindBin(0.01)) + eta->GetBinContent(eta->FindBin(-0.01))) / 2);
2843   Printf("dndeta in |eta| < 0.5 is %f", eta->Integral(eta->FindBin(-0.49), eta->FindBin(0.49)) / (eta->FindBin(0.49) - eta->FindBin(-0.49) + 1));
2844   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));
2845   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));
2846
2847   stats = (TH2*) gFile->Get("fEventStats");
2848   proj = stats->ProjectionX();
2849   gROOT->ProcessLine(".L PrintHist.C");
2850   PrintHist2D(stats);
2851   PrintHist(proj);
2852
2853   Printf("+++ TRIGGER EFFICIENCIES +++");
2854
2855   Printf("INEL = %.1f", 100. * (proj->GetBinContent(1) - stats->GetBinContent(1, 1) - stats->GetBinContent(1, 3)) / proj->GetBinContent(1));
2856   Printf("NSD  = %.1f", 100. * (proj->GetBinContent(2) - stats->GetBinContent(2, 1) - stats->GetBinContent(2, 3)) / proj->GetBinContent(2));
2857   Printf("ND  = %.1f",  100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1) - stats->GetBinContent(3, 3)) / proj->GetBinContent(3));
2858   Printf("SD  = %.1f",  100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1) - stats->GetBinContent(4, 3)) / proj->GetBinContent(4));
2859   Printf("DD  = %.1f",  100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1) - stats->GetBinContent(5, 3)) / proj->GetBinContent(5));
2860
2861   Printf("+++ TRIGGER + VERTEX EFFICIENCIES +++");
2862
2863   Printf("INEL = %.1f", 100. * stats->GetBinContent(1, 4) / proj->GetBinContent(1));
2864   Printf("NSD  = %.1f", 100. * stats->GetBinContent(2, 4) / proj->GetBinContent(2));
2865   Printf("ND  = %.1f",  100. * stats->GetBinContent(3, 4) / proj->GetBinContent(3));
2866   Printf("SD  = %.1f",  100. * stats->GetBinContent(4, 4) / proj->GetBinContent(4));
2867   Printf("DD  = %.1f",  100. * stats->GetBinContent(5, 4) / proj->GetBinContent(5));
2868
2869   for (Int_t i=7; i<=proj->GetNbinsX(); i++)
2870     if (proj->GetBinContent(i) > 0)
2871       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));
2872
2873   //eta->Draw();
2874 }
2875
2876 void TestAsymmetry()
2877 {
2878   loadlibs();
2879
2880   TFile* file2 = TFile::Open("analysis_mc.root");
2881
2882   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
2883   fdNdEtaAnalysis->LoadHistograms();
2884   fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
2885
2886   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy();
2887
2888   hist = (TH1*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
2889   hist2 = (TH1*) hist->Clone("hist2");
2890   for (Int_t x=1; x<=hist->GetNbinsX(); x++)
2891     for (Int_t y=1; y<=hist->GetNbinsY(); y++)
2892       for (Int_t z=1; z<=hist->GetNbinsZ(); z++)
2893       {
2894         Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
2895         hist->SetBinContent(x, y, z, hist2->GetBinContent(hist->GetNbinsX() / 2, TMath::Min(y, hist->GetNbinsY() + 1 - y), z));
2896       }
2897
2898   hist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
2899   for (Int_t x=1; x<=hist->GetNbinsX(); x++)
2900     for (Int_t y=1; y<=hist->GetNbinsY(); y++)
2901       {
2902         //Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
2903         hist->SetBinContent(x, y, hist->GetBinContent(hist->GetNbinsX() / 2, y));
2904       }
2905
2906   fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
2907   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerColor(2);
2908   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetLineColor(2);
2909   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerStyle(5);
2910   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy("SAMEP");
2911 }
2912
2913 void DeltaPhiFromPt(Float_t smearing = 0.005)
2914 {
2915   loadlibs();
2916
2917   TFile::Open("analysis_mc.root");
2918   hist = (TH1*) gFile->Get("dndeta_check_pt");
2919
2920   dPhiHist = new TH1F("dPhiHist", ";#Delta phi", 400, -0.1, 0.1);
2921   dPhiHist2 = new TH1F("dPhiHist2", ";#Delta phi", 400, -0.1, 0.1);
2922
2923   for (Int_t i=1; i<=hist->GetNbinsX(); i++)
2924   {
2925     Float_t pt = hist->GetBinCenter(i);
2926     Float_t deltaPhi = (0.076 - 0.039) / 2 / (pt / 0.15);
2927
2928     if (smearing > 0)
2929     {
2930       gaus = new TF1("mygaus", "gaus(0)", -0.1, 0.1);
2931       gaus->SetParameters(1, -deltaPhi, smearing);
2932
2933       dPhiHist->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
2934
2935       dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
2936       gaus->SetParameters(1, deltaPhi, smearing);
2937       dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
2938     }
2939     else
2940 {
2941 dPhiHist->Fill(deltaPhi, hist->GetBinContent(i) / 2);
2942 dPhiHist2->Fill(deltaPhi, hist->GetBinContent(i) / 2);
2943 dPhiHist2->Fill(-deltaPhi, hist->GetBinContent(i) / 2);
2944 }
2945   }
2946
2947   new TCanvas;
2948   dPhiHist->Draw();
2949   dPhiHist2->SetLineColor(2);
2950   dPhiHist2->Draw("SAME");
2951   gPad->SetLogy();
2952
2953   TFile::Open("trackletsDePhi.root");
2954   //TFile::Open("tmp/correction_maponly-positive.root");
2955   //TFile::Open("tmp/correction_map.root");
2956   //tracklets = (TH1*) gFile->Get(Form("fDeltaPhi_%d", 0));
2957   tracklets = (TH1*) gFile->Get("DePhiPPTracklets");
2958   tracklets->Scale(1.0 / tracklets->GetMaximum() * dPhiHist->GetMaximum());
2959   tracklets->SetLineColor(4);
2960   tracklets->Draw("SAME");
2961 }
2962
2963 void VertexDistributions()
2964 {
2965   loadlibs();
2966
2967   TFile::Open("correction_map.root");
2968   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
2969   if (!dNdEtaCorrection->LoadHistograms())
2970     return;
2971
2972   all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all");
2973   trigger = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("trigger");
2974   vtx = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("vtx");
2975
2976   nottriggered = (TH1*) all->Clone("nottriggered");
2977   nottriggered->Add(trigger, -1);
2978
2979   novertex = (TH1*) trigger->Clone("novertex");
2980   novertex->Add(vtx, -1);
2981
2982   temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
2983   highmult = temphist->ProjectionX("highmult", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
2984   //all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
2985
2986   for (Int_t i=1; i<=trigger->GetNbinsX(); i++)
2987   {
2988     all->SetBinContent(i, all->GetBinContent(i) / all->GetBinWidth(i));
2989     trigger->SetBinContent(i, trigger->GetBinContent(i) / trigger->GetBinWidth(i));
2990     vtx->SetBinContent(i, vtx->GetBinContent(i) / vtx->GetBinWidth(i));
2991     nottriggered->SetBinContent(i, nottriggered->GetBinContent(i) / nottriggered->GetBinWidth(i));
2992     novertex->SetBinContent(i, novertex->GetBinContent(i) / novertex->GetBinWidth(i));
2993     highmult->SetBinContent(i, highmult->GetBinContent(i) / highmult->GetBinWidth(i));
2994   }
2995
2996   new TCanvas;
2997   vtx->SetTitle("");
2998   vtx->SetStats(0);
2999   vtx->DrawCopy("HIST");
3000
3001   all->Scale(1.0 / all->Integral());
3002   nottriggered->Scale(1.0 / nottriggered->Integral());
3003   novertex->Scale(1.0 / novertex->Integral());
3004   highmult->Scale(1.0 / highmult->Integral());
3005
3006   new TCanvas;
3007   all->Draw("HIST");
3008   novertex->SetLineColor(2);
3009   novertex->Draw("HISTSAME");
3010   highmult->SetLineColor(3);
3011   highmult->Draw("HISTSAME");
3012
3013   legend = new TLegend(0.5, 0.5, 0.8, 0.8);
3014   legend->SetFillColor(0);
3015   legend->AddEntry(all, "all");
3016   legend->AddEntry(novertex, "no vertex");
3017   legend->AddEntry(highmult, "mult > 10");
3018   legend->Draw();
3019
3020   new TCanvas;
3021   trigger->Scale(1.0 / trigger->Integral());
3022   vtx->Scale(1.0 / vtx->Integral());
3023
3024   trigger->Divide(vtx);
3025
3026   trigger->Draw();
3027   //vtx->SetLineColor(2);
3028   //vtx->Draw("SAME");
3029
3030   //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
3031   temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram();
3032   //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram();
3033
3034   temphist = (TH2*) gFile->Get("fTemp1");
3035
3036   new TCanvas;
3037   legend = new TLegend(0.7, 0.7, 0.99, 0.99);
3038   legend->SetFillColor(0);
3039
3040   Bool_t first = kTRUE;
3041   for (Int_t i=0; i<20; i+=5)
3042   {
3043     highmult = temphist->ProjectionX("highmult", i+1, i+1+4);
3044     highmult->Rebin(10);
3045     Printf("%f", highmult->Integral());
3046     if (highmult->Integral() <= 0)
3047       continue;
3048
3049     for (Int_t j=1; j<=trigger->GetNbinsX(); j++)
3050       highmult->SetBinContent(j, highmult->GetBinContent(j) / highmult->GetBinWidth(j));
3051
3052     highmult->Scale(1.0 / highmult->Integral());
3053     highmult->SetLineColor((i/5)+1);
3054     highmult->GetYaxis()->SetRangeUser(0, 0.15);
3055     if (first)
3056     {
3057       highmult->DrawCopy();
3058       first = kFALSE;
3059     }
3060     else
3061       highmult->DrawCopy("SAME");
3062     legend->AddEntry(highmult->Clone(), Form("%d <= N <= %d", i, i+4));
3063   }
3064   legend->Draw();
3065
3066 }
3067
3068 void PlotPt1DCorrection()
3069 {
3070   const char* files[] = { "field.root", "field_onlyprim.root", "nofield.root", "nofield_onlyprim.root" };
3071   const char* names[] = { "B: all", "B: primaries", "No B: all", "No B: primaries" };
3072   Int_t colors[] = { 1, 2, 3, 4 };
3073
3074   loadlibs();
3075
3076   dummy = new TH2F("dummy", ";p_{T};correction", 100, 0, 1.4, 100, 0.5, 3);
3077   dummy->SetStats(0);
3078   //dummy->GetYaxis()->SetTitleOffset(1.3);
3079   dummy->Draw();
3080
3081   legend = new TLegend(0.48, 0.57, 0.88, 0.88);
3082   legend->SetFillColor(0);
3083
3084   for (Int_t i=0; i<4; i++)
3085   {
3086     TFile::Open(files[i]);
3087     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
3088     if (!dNdEtaCorrection->LoadHistograms())
3089       return;
3090
3091     hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get1DCorrectionHistogram("z", -9.9, 9.9, -0.79, 0.79);
3092     hist->SetLineColor(colors[i]);
3093     hist->SetLineWidth(2);
3094     hist->SetMarkerColor(colors[i]);
3095     hist->Draw("SAME");
3096
3097     legend->AddEntry(hist, names[i], "L");
3098   }
3099
3100   legend->Draw();
3101 }