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   }
 150 }
 151
 152 void PrintAllInfos()
 153 {
 154   PrintInfo();
 155
 156   Printf("RAW ESD");
 157   TFile::Open("analysis_esd_raw.root");
 158   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 159   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
 160   fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
 161
 162   const Int_t num = 3;
 163   const char* results[] = { "dndeta", "dndetaTr", "dndetaTrVtx" };
 164
 165   TFile::Open("analysis_esd.root");
 166   for (Int_t i=0; i<num; i++)
 167   {
 168     Printf("CORRECTED %s", results[i]);
 169     dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 170     fdNdEtaAnalysis->LoadHistograms(results[i]);
 171     fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
 172   }
 173 }
 174
 175 void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
 176 {
 177   gSystem->Load("libPWG0base");
 178
 179   TFile::Open(esdFile);
 180   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 181   fdNdEtaAnalysisESD->LoadHistograms();
 182
 183   TFile::Open(mcFile);
 184   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
 185   fdNdEtaAnalysisMC->LoadHistograms();
 186   //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 187
 188   for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
 189     fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
 190
 191   fdNdEtaAnalysisESD->DrawHistograms();
 192 }
 193
 194 void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 195 {
 196   gSystem->Load("libPWG0base");
 197
 198   const char* ESDfolder = 0;
 199
 200   if (plot == 0) // all
 201     ESDfolder = "dndeta";
 202   else if (plot == 1) // mb
 203     ESDfolder = "dndeta_mb";
 204   else if (plot == 2) // mb vtx
 205     ESDfolder = "dndeta_mbvtx";
 206
 207   TFile::Open("analysis_esd.root");
 208   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
 209   fdNdEtaAnalysisESD->LoadHistograms();
 210
 211   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 212   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 213
 214   TH2F* hist = 0;
 215
 216   if (plot == 0) // all
 217     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
 218   else if (plot == 1) // mb
 219     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
 220   else if (plot == 2) // mb vtx
 221     hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
 222
 223   TH1* proj = hist->ProjectionX();
 224
 225   TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
 226   for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
 227   {
 228     Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
 229     if (value != 0)
 230     {
 231       printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
 232       vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
 233     }
 234   }
 235
 236   new TCanvas;
 237   vertex->DrawCopy();
 238 }
 239
 240 void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
 241 {
 242   loadlibs();
 243
 244   TFile::Open("analysis_esd.root");
 245   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
 246   fdNdEtaAnalysisESD->LoadHistograms();
 247
 248   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 249   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
 250
 251   //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
 252   //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
 253
 254   TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
 255   TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 256
 257   TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
 258
 259   new TCanvas;
 260   histESD->Draw();
 261
 262   new TCanvas;
 263   histCorr->Draw();
 264
 265   for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
 266     for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
 267       for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
 268       {
 269         Float_t value1 = histESD->GetBinContent(x, y, z);
 270         Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
 271
 272         if (value2 > 0 && value1 > 0)
 273         {
 274           printf("%f %f %f\n", value1, value2, value1 / value2);
 275           diff->Fill(value1 / value2);
 276         }
 277       }
 278
 279   new TCanvas;
 280   diff->Draw();
 281 }
 282
 283 Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
 284 {
 285   Double_t avgMC = 0;
 286   Double_t avgESD = 0;
 287   for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
 288   {
 289     avgMC += histMC->GetBinContent(bin);
 290     avgESD += histESD->GetBinContent(bin);
 291   }
 292   Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
 293
 294   avgMC /= nBins;
 295   avgESD /= nBins;
 296
 297   // deviation when integrate in |eta| < 0.8 between mc and esd
 298   Double_t diffFullRange = (avgMC - avgESD) / avgMC;
 299
 300   Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
 301
 302   return diffFullRange;
 303 }
 304
 305 void dNdEtaNoResolution()
 306 {
 307   loadlibs();
 308
 309   TFile::Open("correction_map.root");
 310
 311   const char* correctionMapFolder = "dndeta_correction";
 312   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
 313   dNdEtaCorrection->LoadHistograms();
 314   dNdEtaCorrection->GetTrack2ParticleCorrection()->PrintInfo(0.3);
 315
 316   TFile::Open("analysis_mc.root");
 317   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
 318   fdNdEtaAnalysis->LoadHistograms();
 319   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, "ESD (no resolution effect) -> MB with vertex");
 320   fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0)->SetMarkerStyle(21);
 321
 322   TFile::Open("analysis_mc.root");
 323   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
 324   fdNdEtaAnalysisMC->LoadHistograms();
 325   fdNdEtaAnalysisMC->Finish(0, 0, AlidNdEtaCorrection::kNone, "MC: MB with vertex");
 326
 327   DrawdNdEtaRatio(fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0), fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(0), "MB with vertex (no resolution effect)", 3);
 328 }
 329
 330 TH1* GetMCHist(const char* folder, Float_t ptCut, const char* tag)
 331 {
 332   loadlibs();
 333
 334   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(folder, folder);
 335   fdNdEtaAnalysis->LoadHistograms();
 336   fdNdEtaAnalysis->Finish(0, ptCut, AlidNdEtaCorrection::kNone, tag);
 337   return fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 338 }
 339
 340 void dNdEtaFinal(Bool_t spd = kTRUE)
 341 {
 342   TFile* file = TFile::Open("analysis_esd.root");
 343   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
 344   TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
 345   Prepare1DPlot(histESD);
 346   Prepare1DPlot(histESDnsd);
 347
 348   TCanvas* canvas = new TCanvas("dNdEtaFinal", "dNdEtaFinal", 600, 600);
 349   gPad->SetTopMargin(0.05);
 350   gPad->SetRightMargin(0.05);
 351   gPad->SetLeftMargin(0.12);
 352   gPad->SetBottomMargin(0.12);
 353   gPad->SetGridx();
 354   gPad->SetGridy();
 355
 356   Float_t etaMax = 1.9;
 357   Float_t histMax = 1.39;
 358   Float_t systErrorValue = 0.023;
 359   Float_t systErrorNSDValue = 0.081;
 360   if (!spd)
 361   {
 362     //etaMax = 1.5;
 363     histMax = 0.99;
 364     systErrorValue = 0.043;
 365     systErrorNSDValue = 0.088;
 366   }
 367
 368   dummy = new TH2F("dummy", ";#eta;dN_{ch}/d#eta", 100, -etaMax, etaMax, 100, 3, 8);
 369   dummy->SetStats(0);
 370   dummy->GetYaxis()->SetTitleOffset(1.3);
 371
 372   histESD->SetMarkerStyle(20);
 373   histESDnsd->SetMarkerStyle(21);
 374   histESDnsd->SetMarkerColor(4);
 375   histESDnsd->SetLineColor(4);
 376   histESD->SetMarkerSize(1.5);
 377   histESDnsd->SetMarkerSize(1.5);
 378
 379   histESD->GetXaxis()->SetRangeUser(-histMax, histMax);
 380   histESDnsd->GetXaxis()->SetRangeUser(-histMax, histMax);
 381
 382   legend = new TLegend(0.3, 0.2, 0.78, 0.4);
 383   legend->SetFillColor(0);
 384   legend->SetTextSize(0.04);
 385   legend->AddEntry(histESD, "Inelastic events", "P");
 386   legend->AddEntry(histESDnsd, "NSD events", "P");
 387
 388   dummy->Draw();
 389
 390   // syst errors.
 391   TH1* systError = (TH1*) histESD->Clone("systError");
 392   for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
 393     systError->SetBinError(i, systError->GetBinContent(i) * systErrorValue);
 394   // change error drawing style
 395   systError->SetFillColor(15);
 396   systError->DrawCopy("SAME E2 ][");
 397
 398   // syst error NSD
 399   for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
 400   {
 401     systError->SetBinContent(i, histESDnsd->GetBinContent(i));
 402     systError->SetBinError(i, systError->GetBinContent(i) * systErrorNSDValue);
 403   }
 404   systError->DrawCopy("SAME E2 ][");
 405
 406   histESD->Draw("SAME");
 407   histESDnsd->Draw("SAME");
 408   legend->Draw();
 409
 410   canvas->SaveAs(Form("%s_dndeta_final.eps", (spd) ? "spd" : "tpc"));
 411 }
 412
 413 void dNdEtaPythiaPhojet()
 414 {
 415   // evtl. deactivate acceptance maps in dNdEtaAnalysis.cxx
 416
 417   loadlibs();
 418
 419   TH1* hist[4];
 420
 421   TFile::Open("LHC08c11_10TeV_0.5T/mb1/spd/analysis_mc.root");
 422   hist[0] =         (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
 423   hist[1] =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
 424
 425   TFile::Open("LHC08c15_10TeV_0.5T_Phojet/mb1/spd/analysis_mc.root");
 426   hist[2] =         (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMCPhojet");
 427   hist[3] =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsdPhojet");
 428
 429   file = TFile::Open("pythia_phojet_dndeta.root", "RECREATE");
 430   for (Int_t i=0; i<4; i++)
 431     hist[i]->Write();
 432   file->Close();
 433 }
 434
 435 void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 436 {
 437   loadlibs();
 438
 439   TFile* file = TFile::Open("analysis_esd.root");
 440
 441   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 442   fdNdEtaAnalysis->LoadHistograms("dndeta");
 443
 444   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
 445   TH1* histESD1 = (TH1*) file->Get("dndeta/dNdEta_corrected_1");
 446   TH1* histESD2 = (TH1*) file->Get("dndeta/dNdEta_corrected_2");
 447   TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
 448   TH1* histESDnsdNoPt = (TH1*) file->Get("dndetaNSD/dNdEta");
 449   TH1* histESDonePart = (TH1*) file->Get("dndetaOnePart/dNdEta_corrected");
 450   if (!histESDonePart)
 451     histESDonePart = new TH1F;
 452   TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
 453   TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
 454   TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
 455   TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
 456   TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
 457   TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
 458
 459   Prepare1DPlot(histESD);
 460   Prepare1DPlot(histESD1);
 461   Prepare1DPlot(histESD2);
 462   Prepare1DPlot(histESDnsd);
 463   Prepare1DPlot(histESDonePart);
 464   Prepare1DPlot(histESDMB);
 465   Prepare1DPlot(histESDMBVtx);
 466
 467   Prepare1DPlot(histESDNoPt);
 468   Prepare1DPlot(histESDMBNoPt);
 469   Prepare1DPlot(histESDMBVtxNoPt);
 470   Prepare1DPlot(histESDMBTracksNoPt);
 471
 472   histESD->SetLineWidth(0);
 473   histESDnsd->SetLineWidth(0);
 474   histESDonePart->SetLineWidth(0);
 475   histESDMB->SetLineWidth(0);
 476   histESDMBVtx->SetLineWidth(0);
 477
 478   histESDNoPt->SetLineWidth(0);
 479   histESDMBNoPt->SetLineWidth(0);
 480   histESDMBVtxNoPt->SetLineWidth(0);
 481
 482   histESD->SetMarkerColor(1);
 483   histESDnsd->SetMarkerColor(6);
 484   histESDonePart->SetMarkerColor(3);
 485   histESDMB->SetMarkerColor(2);
 486   histESDMBVtx->SetMarkerColor(4);
 487
 488   histESD->SetLineColor(1);
 489   histESDnsd->SetLineColor(6);
 490   histESDonePart->SetLineColor(3);
 491   histESDMB->SetLineColor(2);
 492   histESDMBVtx->SetLineColor(4);
 493
 494   histESDNoPt->SetMarkerColor(1);
 495   histESDMBNoPt->SetMarkerColor(2);
 496   histESDMBVtxNoPt->SetMarkerColor(4);
 497   histESDMBTracksNoPt->SetMarkerColor(3);
 498
 499   histESD->SetMarkerStyle(20);
 500   histESDnsd->SetMarkerStyle(29);
 501   histESDonePart->SetMarkerStyle(24);
 502   histESDMB->SetMarkerStyle(21);
 503   histESDMBVtx->SetMarkerStyle(22);
 504
 505   histESDNoPt->SetMarkerStyle(20);
 506   histESDMBNoPt->SetMarkerStyle(21);
 507   histESDMBVtxNoPt->SetMarkerStyle(22);
 508   histESDMBTracksNoPt->SetMarkerStyle(23);
 509
 510   Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 0.89 : 1.99;
 511   Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 1.2 : 2.3;
 512   //Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 0.89 : 1.39;
 513   //Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 1.2 : 1.9;
 514
 515   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 516   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 517   histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 518   histESDnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 519   histESDonePart->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 520
 521   histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 522   histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 523   histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 524   histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 525
 526   Float_t max = TMath::Max(histESDMBVtx->GetMaximum(), histESDMB->GetMaximum());
 527   max = TMath::Max(max, histESD->GetMaximum());
 528   max = TMath::Max(max, histESDMBTracksNoPt->GetMaximum());
 529
 530   TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.4);
 531   legend->SetFillColor(0);
 532   legend->AddEntry(histESDMBVtx, "Triggered, vertex");
 533   legend->AddEntry(histESDMB, "Triggered");
 534   legend->AddEntry(histESD, "All INEL events");
 535   legend->AddEntry(histESDnsd, "All NSD events");
 536   legend->AddEntry(histESDonePart, "One Particle");
 537
 538   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 0, max * 1.1);
 539   dummy->GetYaxis()->SetRangeUser(2.1, max * 1.1);
 540   Prepare1DPlot(dummy);
 541   dummy->SetStats(kFALSE);
 542   dummy->SetXTitle("#eta");
 543   dummy->SetYTitle("dN_{ch}/d#eta");
 544   dummy->GetYaxis()->SetTitleOffset(1);
 545
 546   TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
 547
 548   dummy->DrawCopy();
 549   histESDMBVtx->Draw("SAME");
 550   histESDMB->Draw("SAME");
 551   histESD->Draw("SAME");
 552   histESDnsd->Draw("SAME");
 553   histESDonePart->Draw("SAME");
 554   legend->Draw();
 555
 556   if (save)
 557   {
 558     canvas->SaveAs("dNdEta1.png");
 559     //canvas->SaveAs("dNdEta1.eps");
 560   }
 561
 562   histESD->Fit("pol0", "0", "", -0.49, 0.49);
 563   histESDnsd->Fit("pol0", "0", "", -0.49, 0.49);
 564   histESDonePart->Fit("pol0", "0", "", -0.49, 0.49);
 565   histESDonePart->Fit("pol0", "0", "", -0.99, 0.99);
 566
 567   canvas = new TCanvas("dNdEta1_mirrored", "dNdEta1_mirrored", 500, 500);
 568   canvas->SetGridx();
 569   canvas->SetGridy();
 570
 571   dummy->DrawCopy()->GetXaxis()->SetRangeUser(0, 100);
 572   histESD->DrawCopy("SAME")->SetMarkerStyle(24);
 573   histESDnsd->DrawCopy("SAME")->SetMarkerStyle(24);
 574
 575   graph = new TGraphErrors(histESD);
 576   for (Int_t i=0; i<graph->GetN(); i++)
 577     graph->GetX()[i] *= -1;
 578   graph->SetMarkerStyle(5);
 579   graph->Draw("P SAME");
 580
 581   graph = new TGraphErrors(histESDnsd);
 582   for (Int_t i=0; i<graph->GetN(); i++)
 583     graph->GetX()[i] *= -1;
 584   graph->SetMarkerStyle(5);
 585   graph->SetMarkerColor(histESDnsd->GetMarkerColor());
 586   graph->Draw("P SAME");
 587
 588   canvas = new TCanvas("dNdEta1_ratio", "dNdEta1_ratio", 500, 500);
 589   canvas->SetGridx();
 590   canvas->SetGridy();
 591
 592   dummy_clone = dummy->DrawCopy();
 593   dummy_clone->GetXaxis()->SetRangeUser(0, 100);
 594   dummy_clone->GetYaxis()->SetRangeUser(0.5, 1.5);
 595
 596   graph = new TGraphErrors(histESD);
 597   for (Int_t i=0; i<graph->GetN(); i++)
 598   {
 599     Int_t bin = histESD->GetXaxis()->FindBin(-graph->GetX()[i]);
 600     if (histESD->GetBinContent(bin) > 0 && graph->GetY()[i] > 0)
 601     {
 602       graph->GetEY()[i] = sqrt(graph->GetEY()[i] * graph->GetEY()[i] / graph->GetY()[i] / graph->GetY()[i]
 603         + histESD->GetBinError(bin) * histESD->GetBinError(bin) / histESD->GetBinContent(bin) / histESD->GetBinContent(bin));
 604       graph->GetY()[i] /= histESD->GetBinContent(bin);
 605       graph->GetEY()[i] *= graph->GetY()[i];
 606     }
 607     else
 608       graph->GetY()[i] = 0;
 609   }
 610   graph->SetMarkerStyle(5);
 611   graph->Draw("P SAME");
 612
 613   graph = new TGraphErrors(histESDnsd);
 614   for (Int_t i=0; i<graph->GetN(); i++)
 615   {
 616     Int_t bin = histESDnsd->GetXaxis()->FindBin(-graph->GetX()[i]);
 617     if (histESDnsd->GetBinContent(bin) > 0 && graph->GetY()[i] > 0)
 618     {
 619       graph->GetEY()[i] = sqrt(graph->GetEY()[i] * graph->GetEY()[i] / graph->GetY()[i] / graph->GetY()[i]
 620         + histESDnsd->GetBinError(bin) * histESDnsd->GetBinError(bin) / histESDnsd->GetBinContent(bin) / histESDnsd->GetBinContent(bin));
 621       graph->GetY()[i] /= histESDnsd->GetBinContent(bin);
 622       graph->GetEY()[i] *= graph->GetY()[i];
 623       graph->GetY()[i] += 0.2;
 624     }
 625   }
 626   graph->SetMarkerStyle(5);
 627   graph->SetMarkerColor(histESDnsd->GetMarkerColor());
 628   graph->Draw("P SAME");
 629
 630   canvas = new TCanvas("dNdEta1_vertex", "dNdEta1_vertex", 500, 500);
 631   dummy->DrawCopy();
 632   histESD->DrawCopy("SAME");
 633   histESD1->SetLineColor(2);
 634   histESD1->DrawCopy("SAME");
 635   histESD2->SetLineColor(4);
 636   histESD2->DrawCopy("SAME");
 637
 638   PrintIntegratedDeviation(histESDnsd, histESDMB, "factor MB / NSD");
 639
 640   if (onlyESD)
 641     return;
 642
 643   loadlibs();
 644
 645   TFile* file2 = TFile::Open("analysis_mc.root");
 646
 647   TH1* histMCTrVtx =       (TH1*) GetMCHist("dndetaTrVtx", -1, "MC: MB with vertex")->Clone("histMCTrVtx");
 648   TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
 649
 650   TH1* histMC =            (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
 651   TH1* histMCTr =          (TH1*) GetMCHist("dndetaTr", -1, "MC: minimum bias")->Clone("histMCTr");
 652   TH1* histMCnsd =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
 653   TH1* histMConePart =     (TH1*) GetMCHist("dndetaOnePart", -1, "MC: OnePart")->Clone("histMConePart");
 654
 655   TH1* histMCPtCut =       (TH1*) GetMCHist("dndeta", 0.151, "MC: full inelastic, pt cut")->Clone("histMCPtCut");
 656   TH1* histMCTrPtCut =     (TH1*) GetMCHist("dndetaTr", 0.151, "MC: minimum bias, pt cut")->Clone("histMCTrPtCut");
 657   TH1* histMCTrVtxPtCut =  (TH1*) GetMCHist("dndetaTrVtx", 0.151, "MC: MB with vertex, pt cut")->Clone("histMCTrVtxPtCut");
 658   TH1* histMCnsdNoPt =     (TH1*) GetMCHist("dndetaNSD", 0.151, "MC: NSD, put cut")->Clone("histMCnsdNoPt");
 659   TH1* histMCTracksPtCut = (TH1*) GetMCHist("dndetaTracks", 0.151, "MC: Tracks w/o resolution effect, pt cut")->Clone("histMCTracksPtCut");
 660
 661   Prepare1DPlot(histMC);
 662   Prepare1DPlot(histMCnsd);
 663   Prepare1DPlot(histMCTr);
 664   Prepare1DPlot(histMCTrVtx);
 665
 666   Prepare1DPlot(histMCPtCut);
 667   Prepare1DPlot(histMCTrPtCut);
 668   Prepare1DPlot(histMCTrVtxPtCut);
 669   Prepare1DPlot(histMCTracksPtCut);
 670
 671   histMC->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 672   histMCnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 673   histMCTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 674   histMCTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 675
 676   histMCPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 677   histMCTrPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 678   histMCTrVtxPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 679   histMCTracksPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 680
 681   histMC->SetLineColor(1);
 682   histMCnsd->SetLineColor(6);
 683   histMCTr->SetLineColor(2);
 684   histMCTrVtx->SetLineColor(4);
 685
 686   histMCPtCut->SetLineColor(1);
 687   histMCTrPtCut->SetLineColor(2);
 688   histMCTrVtxPtCut->SetLineColor(4);
 689   if (histMCTracksPtCut)
 690     histMCTracksPtCut->SetLineColor(3);
 691
 692   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
 693
 694   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
 695   dummy2->GetYaxis()->SetRangeUser(0, max * 1.1);
 696
 697   dummy2->DrawCopy();
 698   histMC->Draw("SAME");
 699   histMCnsd->Draw("SAME");
 700   histMCTr->Draw("SAME");
 701   histMCTrVtx->Draw("SAME");
 702   histESD->Draw("SAME");
 703   histESDnsd->Draw("SAME");
 704   histESDMB->Draw("SAME");
 705   histESDMBVtx->Draw("SAME");
 706   histESDNoPt->Draw("SAME");
 707   histESDMBNoPt->Draw("SAME");
 708   histESDMBVtxNoPt->Draw("SAME");
 709   histESDMBTracksNoPt->Draw("SAME");
 710   histMCPtCut->Draw("SAME");
 711   histMCTrPtCut->Draw("SAME");
 712   histMCTrVtxPtCut->Draw("SAME");
 713   if (histMCTracksPtCut)
 714     histMCTracksPtCut->Draw("SAME");
 715
 716   if (save)
 717   {
 718     canvas2->SaveAs("dNdEta2.gif");
 719     canvas2->SaveAs("dNdEta2.eps");
 720   }
 721
 722   TH1* ratio = (TH1*) DrawdNdEtaRatio(histESD, histMC, "full_inelastic", etaPlotLimit)->Clone();
 723   TH1* ratioTr = (TH1*) DrawdNdEtaRatio(histESDMB, histMCTr, "triggered", etaPlotLimit)->Clone();
 724   TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
 725   TH1* ratioTrVtxNoPt = (TH1*) DrawdNdEtaRatio(histESDMBVtxNoPt, histMCTrVtxPtCut, "triggered_vertex_nopt", etaPlotLimit)->Clone();
 726   TH1* ratioNSD = (TH1*) DrawdNdEtaRatio(histESDnsd, histMCnsd, "NSD", etaPlotLimit)->Clone();
 727   TH1* ratioOnePart = (TH1*) DrawdNdEtaRatio(histESDonePart, histMConePart, "OnePart", etaPlotLimit)->Clone();
 728
 729   // draw ratios of single steps
 730   c7 = new TCanvas("all_ratios", "all_ratios", 600, 600);
 731   c7->SetRightMargin(0.05);
 732   c7->SetTopMargin(0.05);
 733   c7->SetGridx();
 734   c7->SetGridy();
 735
 736   ratioTrVtxNoPt->SetMarkerStyle(20);
 737   ratioTrVtx->SetMarkerStyle(21);
 738   ratioTr->SetMarkerStyle(23);
 739   ratio->SetMarkerStyle(22);
 740   ratioNSD->SetMarkerStyle(26);
 741
 742   ratioTrVtxNoPt->SetMarkerSize(2);
 743   ratioTrVtx->SetMarkerSize(2);
 744   ratioTr->SetMarkerSize(2);
 745   ratio->SetMarkerSize(2);
 746   ratioNSD->SetMarkerSize(2);
 747
 748   ratioTrVtxNoPt->SetMarkerColor(1);
 749   ratioTrVtx->SetMarkerColor(2);
 750   ratioTr->SetMarkerColor(4);
 751   ratio->SetMarkerColor(2);
 752   ratioNSD->SetMarkerColor(1);
 753
 754   ratioTrVtxNoPt->SetLineColor(1);
 755   ratioTrVtx->SetLineColor(2);
 756   ratioTr->SetLineColor(4);
 757   ratio->SetLineColor(2);
 758   ratioNSD->SetLineColor(1);
 759
 760   legend7 = new TLegend(0.13, 0.7, 0.94, 0.9);
 761   legend7->SetFillColor(0);
 762   legend7->SetTextSize(0.035);
 763   legend7->SetNColumns(2);
 764
 765   flat = new TF1("flat", "-1", -5, 5);
 766   ratioTrVtxNoPt->Add(flat);
 767   ratioTrVtx->Add(flat);
 768   ratioTr->Add(flat);
 769   ratio->Add(flat);
 770   ratioNSD->Add(flat);
 771
 772   ratioTrVtxNoPt->Scale(100);
 773   ratioTrVtx->Scale(100);
 774   ratioTr->Scale(100);
 775   ratio->Scale(100);
 776   ratioNSD->Scale(100);
 777
 778   ratio->Add(ratioTr, -1);
 779   ratioNSD->Add(ratioTr, -1);
 780   ratioTr->Add(ratioTrVtx, -1);
 781   ratioTrVtx->Add(ratioTrVtxNoPt, -1);
 782
 783   legend7->AddEntry(ratioTrVtxNoPt, "Track-to-particle", "P");
 784   legend7->AddEntry(ratio, "Trigger-bias INEL", "P");
 785   legend7->AddEntry(ratioTr, "Vertex-reconstruction", "P");
 786   legend7->AddEntry(ratioNSD, "Trigger-bias NSD", "P");
 787   if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
 788     legend7->AddEntry(ratioTrVtx, "p_{T} cut-off", "P");
 789
 790   TH1* dummy7 = new TH2F("dummy7", ";#eta;Deviation in %", 100, -etaPlotLimit, etaPlotLimit, 100, -5, 7);
 791   dummy7->SetStats(0);
 792   dummy7->Draw();
 793
 794   ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 795   ratioTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 796   ratioTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 797   ratioTrVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 798   ratioNSD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 799
 800   ratio->Draw("HIST EP SAME");
 801   ratioTr->Draw("HIST EP SAME");
 802   if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
 803     ratioTrVtx->Draw("HIST EP SAME");
 804   ratioTrVtxNoPt->Draw("HIST EP SAME");
 805   ratioNSD->Draw("HIST EP SAME");
 806   legend7->Draw();
 807
 808   //c7->SaveAs("ratios.eps");
 809
 810   new TCanvas;
 811   dummy2->DrawCopy();
 812   histMCnsd->Draw("SAME");
 813   histESDnsd->Draw("SAME");
 814
 815   ratio = (TH1*) histMC->Clone("ratio");
 816   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
 817
 818   ratio->Divide(histESD);
 819   ratioNoPt->Divide(histESDNoPt);
 820
 821   ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 822
 823   ratio->SetLineColor(1);
 824   ratioNoPt->SetLineColor(2);
 825
 826   Double_t average = 0;       // average deviation from 1 in ratio (depends on the number of bins if statistical)
 827   for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
 828     average += TMath::Abs(ratio->GetBinContent(bin) - 1);
 829   Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
 830   average /= nBins;
 831   Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
 832
 833   PrintIntegratedDeviation(histMC, histESD, "all events");
 834   PrintIntegratedDeviation(histMCnsd, histESDnsd, "all events (NSD)");
 835   PrintIntegratedDeviation(histMConePart, histESDonePart, "all events (INEL>0)");
 836   PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
 837   PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
 838   PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
 839   PrintIntegratedDeviation(histMCnsdNoPt, histESDnsdNoPt, "all events (NSD) (no pt corr)");
 840   PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
 841   PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
 842   PrintIntegratedDeviation(histESD, histESDNoPt, "pt cut off correction");
 843
 844   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 600, 600);
 845   canvas3->Range(0, 0, 1, 1);
 846   //canvas3->Divide(1, 2, 0, 0);
 847
 848   //canvas3->cd(1);
 849   TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
 850   pad1->SetTopMargin(0.05);
 851   pad1->SetLeftMargin(0.13);
 852   pad1->Draw();
 853
 854   TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
 855   pad2->SetLeftMargin(0.13);
 856   pad2->Draw();
 857
 858   pad1->SetRightMargin(0.01);
 859   pad2->SetRightMargin(0.01);
 860
 861   // no border between them
 862   pad1->SetBottomMargin(0);
 863   pad2->SetTopMargin(0);
 864
 865   pad1->cd();
 866   pad1->SetGridx();
 867   pad1->SetGridy();
 868
 869   legend->AddEntry(histMC, "MC prediction");
 870
 871   dummy->GetXaxis()->SetLabelSize(0.08);
 872   dummy->GetYaxis()->SetLabelSize(0.08);
 873   dummy->GetXaxis()->SetTitleSize(0.08);
 874   dummy->GetYaxis()->SetTitleSize(0.08);
 875   dummy->GetYaxis()->SetTitleOffset(0.8);
 876   dummy->DrawCopy();
 877   histESDMBVtx->Draw("SAME");
 878   histESDMB->Draw("SAME");
 879   histESD->Draw("SAME");
 880   histMC->Draw("SAME");
 881
 882   legend->SetTextSize(0.08);
 883   legend->Draw();
 884
 885   pad2->cd();
 886   pad2->SetBottomMargin(0.15);
 887   //pad2->SetGridx();
 888   //pad2->SetGridy();
 889
 890   Float_t minR = 0.91; //TMath::Min(0.961, ratio->GetMinimum() * 0.95);
 891   Float_t maxR = 1.09; //TMath::Max(1.049, ratio->GetMaximum() * 1.05);
 892
 893   TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
 894   dummy3.SetStats(kFALSE);
 895   for (Int_t i=1; i<=100; ++i)
 896     dummy3.SetBinContent(i, 1);
 897   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
 898   dummy3.SetLineWidth(2);
 899   dummy3.GetXaxis()->SetLabelSize(0.08);
 900   dummy3.GetYaxis()->SetLabelSize(0.08);
 901   dummy3.GetXaxis()->SetTitleSize(0.08);
 902   dummy3.GetYaxis()->SetTitleSize(0.08);
 903   dummy3.GetYaxis()->SetTitleOffset(0.8);
 904   dummy3.DrawCopy();
 905
 906   ratio->Draw("SAME");
 907
 908   //pad2->Draw();
 909
 910   canvas3->Modified();
 911
 912   if (save)
 913   {
 914     canvas3->SaveAs("dNdEta.gif");
 915     canvas3->SaveAs("dNdEta.eps");
 916   }
 917
 918   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
 919
 920   ratio->Draw();
 921   ratioNoPt->Draw("SAME");
 922
 923   TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
 924   legend->SetFillColor(0);
 925   legend->AddEntry(ratio, "mc/esd");
 926   legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
 927   legend->Draw();
 928 }
 929
 930 void CompareTwodNdEta(const char* fileName1, const char* fileName2, Bool_t errorsCorrelated = kFALSE)
 931 {
 932   c = new TCanvas;
 933
 934   c->SetGridx();
 935   c->SetGridy();
 936
 937   hist = new TH2F("dummy", ";#eta;dN_{ch}/d#eta", 100, -2.5, 2.5, 100, 0, 8);
 938   hist->SetStats(0);
 939   hist->DrawCopy();//->GetYaxis()->SetRangeUser(2, 4.5);
 940
 941   l = new TLegend(0.2, 0.13, 0.8, 0.35);
 942   l->SetNColumns(2);
 943   l->SetFillColor(0);
 944
 945   TH1* histESD[2];
 946   TH1* histESDnsd[2];
 947
 948   for (Int_t i=0; i<2; i++)
 949   {
 950     if (i == 0)
 951       file = TFile::Open(fileName1);
 952     if (i == 1)
 953     {
 954       if (fileName2 == 0)
 955         break;
 956       file = TFile::Open(fileName2);
 957     }
 958
 959     histESD[i] = (TH1*) file->Get("dndeta/dNdEta_corrected");
 960     histESDnsd[i] = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
 961
 962     histESD[i]->SetMarkerStyle(20 + i*4);
 963     histESDnsd[i]->SetMarkerStyle(21 + i*4);
 964
 965     histESD[i]->SetMarkerColor(i+1);
 966     histESD[i]->SetLineColor(i+1);
 967     histESDnsd[i]->SetMarkerColor(i+1);
 968     histESDnsd[i]->SetLineColor(i+1);
 969
 970     histESD[i]->DrawCopy("SAME");
 971     histESDnsd[i]->DrawCopy("SAME");
 972
 973     l->AddEntry(histESD[i], Form("Data %d INEL", i), "P");
 974     l->AddEntry(histESDnsd[i], Form("Data %d NSD", i), "P");
 975   }
 976
 977   if (0)
 978   {
 979     TGraphErrors *gre = new TGraphErrors(16);
 980     gre->SetFillColor(4);
 981     gre->SetMarkerColor(4);
 982     gre->SetMarkerStyle(26);
 983     gre->SetPoint(0,0.125,3.14);
 984     gre->SetPointError(0,0,0.07);
 985     gre->SetPoint(1,0.375,3.04);
 986     gre->SetPointError(1,0,0.07);
 987     gre->SetPoint(2,0.625,3.17);
 988     gre->SetPointError(2,0,0.07);
 989     gre->SetPoint(3,0.875,3.33);
 990     gre->SetPointError(3,0,0.07);
 991     gre->SetPoint(4,1.125,3.33);
 992     gre->SetPointError(4,0,0.07);
 993     gre->SetPoint(5,1.375,3.53);
 994     gre->SetPointError(5,0,0.07);
 995     gre->SetPoint(6,1.625,3.46);
 996     gre->SetPointError(6,0,0.07);
 997     gre->SetPoint(7,1.875,3.41);
 998     gre->SetPointError(7,0,0.07);
 999     gre->SetPoint(8,-0.125,3.14);
1000     gre->SetPointError(8,0,0.07);
1001     gre->SetPoint(9,-0.375,3.04);
1002     gre->SetPointError(9,0,0.07);
1003     gre->SetPoint(10,-0.625,3.17);
1004     gre->SetPointError(10,0,0.07);
1005     gre->SetPoint(11,-0.875,3.33);
1006     gre->SetPointError(11,0,0.07);
1007     gre->SetPoint(12,-1.125,3.33);
1008     gre->SetPointError(12,0,0.07);
1009     gre->SetPoint(13,-1.375,3.53);
1010     gre->SetPointError(13,0,0.07);
1011     gre->SetPoint(14,-1.625,3.46);
1012     gre->SetPointError(14,0,0.07);
1013     gre->SetPoint(15,-1.875,3.41);
1014     gre->SetPointError(15,0,0.07);
1015     gre->Draw("p");
1016
1017     l->AddEntry(gre, "UA5 INEL", "P");
1018
1019     gre = new TGraphErrors(16);
1020     gre->SetMarkerColor(4);
1021     gre->SetFillColor(4);
1022     gre->SetMarkerStyle(22);
1023     gre->SetPoint(0,0.125,3.48);
1024     gre->SetPointError(0,0,0.07);
1025     gre->SetPoint(1,0.375,3.38);
1026     gre->SetPointError(1,0,0.07);
1027     gre->SetPoint(2,0.625,3.52);
1028     gre->SetPointError(2,0,0.07);
1029     gre->SetPoint(3,0.875,3.68);
1030     gre->SetPointError(3,0,0.07);
1031     gre->SetPoint(4,1.125,3.71);
1032     gre->SetPointError(4,0,0.07);
1033     gre->SetPoint(5,1.375,3.86);
1034     gre->SetPointError(5,0,0.07);
1035     gre->SetPoint(6,1.625,3.76);
1036     gre->SetPointError(6,0,0.07);
1037     gre->SetPoint(7,1.875,3.66);
1038     gre->SetPointError(7,0,0.07);
1039     gre->SetPoint(8,-0.125,3.48);
1040     gre->SetPointError(8,0,0.07);
1041     gre->SetPoint(9,-0.375,3.38);
1042     gre->SetPointError(9,0,0.07);
1043     gre->SetPoint(10,-0.625,3.52);
1044     gre->SetPointError(10,0,0.07);
1045     gre->SetPoint(11,-0.875,3.68);
1046     gre->SetPointError(11,0,0.07);
1047     gre->SetPoint(12,-1.125,3.71);
1048     gre->SetPointError(12,0,0.07);
1049     gre->SetPoint(13,-1.375,3.86);
1050     gre->SetPointError(13,0,0.07);
1051     gre->SetPoint(14,-1.625,3.76);
1052     gre->SetPointError(14,0,0.07);
1053     gre->SetPoint(15,-1.875,3.66);
1054     gre->SetPointError(15,0,0.07);
1055     gre->Draw("p");
1056
1057     l->AddEntry(gre, "UA5 NSD", "P");
1058   }
1059
1060   l->Draw();
1061
1062   if (fileName2 == 0)
1063     return;
1064
1065   new TCanvas;
1066   gPad->SetGridx();
1067   gPad->SetGridy();
1068
1069   if (errorsCorrelated)
1070   {
1071     for (Int_t i=1; i<=histESD[1]->GetNbinsX(); i++)
1072     {
1073       histESD[1]->SetBinError(i, 0);
1074       histESDnsd[1]->SetBinError(i, 0);
1075     }
1076   }
1077
1078   histESD[0]->Divide(histESD[0], histESD[1]);
1079   histESDnsd[0]->Divide(histESDnsd[0], histESDnsd[1]);
1080
1081   for (Int_t i=1; i<=histESD[1]->GetNbinsX(); i++)
1082     histESDnsd[0]->SetBinContent(i, histESDnsd[0]->GetBinContent(i) + 0.2);
1083
1084   hist->DrawCopy()->GetYaxis()->SetRangeUser(0.8, 1.4);
1085   histESD[0]->Draw("SAME");
1086   histESDnsd[0]->Draw("SAME");
1087 }
1088
1089 TH1* DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
1090 {
1091   TCanvas* canvas3 = new TCanvas(name, name, 600, 600);
1092   canvas3->Range(0, 0, 1, 1);
1093
1094   TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
1095   pad1->Draw();
1096
1097   TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
1098   pad2->Draw();
1099
1100   pad1->SetRightMargin(0.01);
1101   pad2->SetRightMargin(0.01);
1102   pad1->SetTopMargin(0.05);
1103   pad1->SetLeftMargin(0.13);
1104   pad2->SetLeftMargin(0.13);
1105   pad2->SetBottomMargin(0.15);
1106
1107   // no border between them
1108   pad1->SetBottomMargin(0);
1109   pad2->SetTopMargin(0);
1110
1111   pad1->cd();
1112   pad1->SetGridx();
1113   pad1->SetGridy();
1114
1115   TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.3);
1116   legend->SetFillColor(0);
1117   legend->AddEntry(corr, "Corrected");
1118   legend->AddEntry(mc, "MC prediction");
1119   legend->SetTextSize(0.08);
1120
1121   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 2.7, corr->GetMaximum() * 1.1);
1122   Prepare1DPlot(dummy);
1123   dummy->SetStats(kFALSE);
1124   dummy->SetXTitle("#eta");
1125   dummy->SetYTitle("dN_{ch}/d#eta");
1126   dummy->GetYaxis()->SetTitleOffset(1);
1127
1128   dummy->GetXaxis()->SetLabelSize(0.08);
1129   dummy->GetYaxis()->SetLabelSize(0.08);
1130   dummy->GetXaxis()->SetTitleSize(0.08);
1131   dummy->GetYaxis()->SetTitleSize(0.08);
1132   dummy->GetYaxis()->SetTitleOffset(0.8);
1133   dummy->DrawCopy();
1134
1135   corr->Draw("SAME");
1136   mc->Draw("SAME");
1137
1138   legend->Draw();
1139
1140   pad2->cd();
1141   pad2->SetBottomMargin(0.15);
1142   //pad2->SetGridx();
1143   //pad2->SetGridy();
1144
1145   TH1* ratio = (TH1*) mc->Clone("ratio");
1146   ratio->Divide(corr);
1147
1148   Float_t minR = TMath::Min(0.91, ratio->GetMinimum() * 0.95);
1149   Float_t maxR = TMath::Max(1.09, ratio->GetMaximum() * 1.05);
1150
1151   TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
1152   dummy3.SetStats(kFALSE);
1153   for (Int_t i=1; i<=100; ++i)
1154         dummy3.SetBinContent(i, 1);
1155   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
1156   dummy3.SetLineWidth(2);
1157   dummy3.GetXaxis()->SetLabelSize(0.08);
1158   dummy3.GetYaxis()->SetLabelSize(0.08);
1159   dummy3.GetXaxis()->SetTitleSize(0.08);
1160   dummy3.GetYaxis()->SetTitleSize(0.08);
1161   dummy3.GetYaxis()->SetTitleOffset(0.8);
1162   dummy3.DrawCopy();
1163
1164   ratio->Draw("SAME");
1165
1166   canvas3->Modified();
1167
1168   return ratio;
1169 }
1170
1171 void ptSpectrum()
1172 {
1173   TFile* file = TFile::Open("analysis_esd.root");
1174   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
1175
1176   TFile* file2 = TFile::Open("analysis_mc.root");
1177   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
1178
1179   TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
1180   InitPad();
1181   gPad->SetLogy();
1182
1183   Prepare1DPlot(histMC);
1184   Prepare1DPlot(histESD);
1185
1186   histESD->SetTitle("");
1187   histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1188   histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
1189
1190   histMC->SetLineColor(kBlue);
1191   histESD->SetLineColor(kRed);
1192
1193   histESD->GetYaxis()->SetTitleOffset(1.5);
1194   histESD->GetXaxis()->SetRangeUser(0, 4.9999);
1195
1196   histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
1197
1198   histESD->Draw();
1199   histMC->Draw("SAME");
1200
1201   canvas->SaveAs("ptSpectrum.gif");
1202   canvas->SaveAs("ptSpectrum.eps");
1203 }
1204
1205 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1206 {
1207   gSystem->Load("libPWG0base");
1208
1209   TFile::Open(fileName);
1210   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1211   dNdEtaCorrection->LoadHistograms();
1212
1213   TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
1214   TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
1215
1216   Prepare2DPlot(corrTrigger);
1217   corrTrigger->SetTitle("b) Trigger bias correction");
1218
1219   Prepare2DPlot(corrVtx);
1220   corrVtx->SetTitle("a) Vertex reconstruction correction");
1221
1222   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
1223   corrVtx->GetYaxis()->SetTitle("Multiplicity");
1224
1225   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
1226   canvas->Divide(2, 1);
1227
1228   canvas->cd(1);
1229   InitPadCOLZ();
1230   corrVtx->DrawCopy("COLZ");
1231
1232   canvas->cd(2);
1233   InitPadCOLZ();
1234   corrTrigger->DrawCopy("COLZ");
1235
1236   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
1237   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
1238
1239   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
1240   canvas->Divide(2, 1);
1241
1242   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
1243   corrVtx->GetYaxis()->SetRangeUser(0, 5);
1244
1245   canvas->cd(1);
1246   InitPadCOLZ();
1247   corrVtx->DrawCopy("COLZ");
1248
1249   canvas->cd(2);
1250   InitPadCOLZ();
1251   corrTrigger->DrawCopy("COLZ");
1252
1253   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
1254   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
1255 }
1256
1257 void TriggerBias(const char* fileName = "correction_map.root")
1258 {
1259   TFile* file = TFile::Open(fileName);
1260
1261   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
1262
1263   Prepare2DPlot(corr);
1264   corr->SetTitle("Trigger bias correction");
1265
1266   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
1267   InitPadCOLZ();
1268   corr->DrawCopy("COLZ");
1269
1270   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
1271   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
1272
1273   corr->GetYaxis()->SetRangeUser(0, 5);
1274
1275   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
1276   InitPadCOLZ();
1277   corr->DrawCopy("COLZ");
1278
1279   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
1280   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
1281 }
1282
1283 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1284 {
1285   gSystem->Load("libPWG0base");
1286
1287   TFile* file = TFile::Open(fileName);
1288   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1289   dNdEtaCorrection->LoadHistograms();
1290
1291   TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->Get1DCorrection("x");
1292   TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->Get1DCorrection("y", -5, 5);
1293
1294   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
1295   canvas->Divide(2, 1);
1296
1297   canvas->cd(1);
1298   InitPad();
1299
1300   Prepare1DPlot(hist);
1301   hist->SetTitle("");
1302   hist->GetYaxis()->SetTitle("correction factor");
1303   hist->GetYaxis()->SetRangeUser(1, 1.5);
1304   hist->GetYaxis()->SetTitleOffset(1.6);
1305   hist->Draw();
1306
1307   canvas->cd(2);
1308   InitPad();
1309
1310   Prepare1DPlot(hist2);
1311   hist2->SetTitle("");
1312   hist2->GetYaxis()->SetTitle("correction factor");
1313   hist2->GetXaxis()->SetRangeUser(0, 5);
1314   hist2->GetYaxis()->SetTitleOffset(1.6);
1315   hist2->GetXaxis()->SetTitle("multiplicity");
1316   hist2->Draw();
1317
1318   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1319   pave->SetFillColor(0);
1320   pave->AddText("|z| < 5 cm");
1321   pave->Draw();
1322
1323   Float_t triggerEff = 100.0 / hist2->GetBinContent(1);
1324   Printf("trigger eff in 0 bin is: %.2f %%", triggerEff);
1325
1326   return;
1327
1328   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1329   //new TCanvas;
1330   //hist2->Draw();
1331
1332   Printf("vertex reco eff in 0 bin is: %.2f %%", 100.0 / hist2->GetBinContent(1));
1333
1334   Printf("combined efficiency is %.2f %%", triggerEff / hist2->GetBinContent(1));
1335 }
1336
1337 void VtxRecon()
1338 {
1339   TFile* file = TFile::Open("correction_map.root");
1340
1341   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
1342
1343   Prepare2DPlot(corr);
1344   corr->SetTitle("Vertex reconstruction correction");
1345
1346   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
1347   InitPadCOLZ();
1348   corr->DrawCopy("COLZ");
1349
1350   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
1351   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
1352
1353   corr->GetYaxis()->SetRangeUser(0, 5);
1354
1355   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
1356   InitPadCOLZ();
1357   corr->DrawCopy("COLZ");
1358
1359   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
1360   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
1361 }
1362
1363 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1364 {
1365   gSystem->Load("libPWG0base");
1366
1367   TFile* file = TFile::Open(fileName);
1368   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1369   dNdEtaCorrection->LoadHistograms();
1370
1371   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
1372   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1373
1374   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
1375   canvas->Divide(2, 1);
1376
1377   canvas->cd(1);
1378   InitPad();
1379
1380   Prepare1DPlot(hist);
1381   hist->SetTitle("");
1382   hist->GetYaxis()->SetTitle("correction factor");
1383   hist->GetYaxis()->SetRangeUser(1, 1.8);
1384   hist->GetYaxis()->SetTitleOffset(1.6);
1385   hist->DrawCopy();
1386
1387   canvas->cd(2);
1388   InitPad();
1389
1390   Prepare1DPlot(hist2);
1391   hist2->SetTitle("");
1392   hist2->GetYaxis()->SetTitle("correction factor");
1393   hist2->GetXaxis()->SetRangeUser(0, 20);
1394   hist2->GetYaxis()->SetTitleOffset(1.6);
1395   hist2->GetXaxis()->SetTitle("multiplicity");
1396   hist2->Draw();
1397
1398   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1399   pave->SetFillColor(0);
1400   pave->AddText("|z| < 10 cm");
1401   pave->Draw();
1402
1403   canvas->SaveAs("VtxRecon1D.eps");
1404
1405   Correction1DCreatePlots(fileName, folderName, 9.9, 2);
1406
1407   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1408   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1409
1410   Prepare1DPlot(corrX);
1411   Prepare1DPlot(corrZ);
1412
1413   corrX->GetYaxis()->SetTitleOffset(1.5);
1414   corrZ->GetYaxis()->SetTitleOffset(1.5);
1415
1416   corrX->SetTitle("a) z projection");
1417   corrZ->SetTitle("b) p_{T} projection");
1418
1419   corrX->GetYaxis()->SetTitle("Correction factor");
1420   corrZ->GetYaxis()->SetTitle("Correction factor");
1421
1422   corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
1423
1424   TString canvasName;
1425   canvasName.Form("VtxRecon1D_Track");
1426   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
1427   canvas->Divide(2, 1);
1428
1429   canvas->cd(1);
1430   InitPad();
1431   corrX->DrawCopy();
1432
1433   canvas->cd(2);
1434   InitPad();
1435   gPad->SetLogx();
1436   corrZ->Draw();
1437
1438   canvas->SaveAs("VtxRecon1D_Track.eps");
1439   canvas->SaveAs("VtxRecon1D_Track.gif");
1440 }
1441
1442 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
1443 {
1444   gSystem->Load("libPWG0base");
1445
1446   TFile::Open(fileName);
1447   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1448   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
1449
1450   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1451   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1452
1453   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1454   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1455
1456   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1457   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1458   gene->GetXaxis()->SetRangeUser(-10, 10);
1459   meas->GetXaxis()->SetRangeUser(-10, 10);
1460
1461   Float_t eff1 = gene->Integral() / meas->Integral();
1462   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1463
1464   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
1465
1466   gene->GetZaxis()->SetRangeUser(0.3, 10);
1467   meas->GetZaxis()->SetRangeUser(0.3, 10);
1468
1469   Float_t eff2 = gene->Integral() / meas->Integral();
1470   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1471
1472   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
1473
1474   gene->GetZaxis()->SetRangeUser(0.3, 1);
1475   meas->GetZaxis()->SetRangeUser(0.3, 1);
1476
1477   Float_t eff3 = gene->Integral() / meas->Integral();
1478   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1479
1480   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
1481 }
1482
1483 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)
1484 {
1485   if (correctionType2 == -1)
1486     correctionType2 = correctionType;
1487
1488   TFile::Open(fileName);
1489   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1490   dNdEtaCorrection->LoadHistograms();
1491
1492   TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
1493   TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType2)->GetTrackCorrection()->GetMeasuredHistogram();
1494
1495   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1496   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1497   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1498   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1499   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kFALSE);
1500   gene->GetYaxis()->SetRange(0, 0);
1501   meas->GetYaxis()->SetRange(0, 0);
1502
1503   gene->GetXaxis()->SetRangeUser(-9.9, 9.9);
1504   meas->GetXaxis()->SetRangeUser(-9.9, 9.9);
1505   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kFALSE);
1506   gene->GetZaxis()->SetRange(0, 0);
1507   meas->GetZaxis()->SetRange(0, 0);
1508
1509   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1510   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1511   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kFALSE);
1512 }
1513
1514 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)
1515 {
1516   gSystem->Load("libPWG0base");
1517
1518   Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType, correctionType2);
1519
1520   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1521   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1522   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1523
1524   Prepare1DPlot(corrX);
1525   Prepare1DPlot(corrY);
1526   Prepare1DPlot(corrZ);
1527
1528   /*
1529   corrX->SetTitle("a) z projection");
1530   corrY->SetTitle("b) #eta projection");
1531   corrZ->SetTitle("c) p_{T} projection");
1532   */
1533
1534   corrX->SetTitle("");
1535   corrY->SetTitle("");
1536   corrZ->SetTitle("");
1537
1538   corrX->SetTitleSize(0.06, "xyz");
1539   corrX->SetLabelSize(0.06, "xyz");
1540   corrY->SetTitleSize(0.06, "xyz");
1541   corrY->SetLabelSize(0.06, "xyz");
1542   corrZ->SetTitleSize(0.06, "xyz");
1543   corrZ->SetLabelSize(0.06, "xyz");
1544
1545   corrX->GetYaxis()->SetTitle("Correction factor");
1546   corrY->GetYaxis()->SetTitle("Correction factor");
1547   corrZ->GetYaxis()->SetTitle("Correction factor");
1548   //corrX->GetYaxis()->SetTitleOffset(1.7);
1549   //corrY->GetYaxis()->SetTitleOffset(1.7);
1550   //corrZ->GetYaxis()->SetTitleOffset(1.7);
1551   corrX->GetYaxis()->SetRangeUser(0.8, 1.5);
1552   corrY->GetYaxis()->SetRangeUser(0.8, 1.5);
1553   corrZ->GetYaxis()->SetRangeUser(0.8, 1.5);
1554
1555   corrZ->GetXaxis()->SetRangeUser(0.11, upperPtLimit);
1556
1557   TString canvasName;
1558   canvasName.Form(Form("Correction1D_%d_%s_%f", correctionType, fileName, upperPtLimit));
1559   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1560   canvas->Divide(3, 1);
1561
1562   TLatex* Tl = new TLatex;
1563   Tl->SetTextSize(0.06);
1564   Tl->SetBit(TLatex::kTextNDC);
1565
1566   canvas->cd(1);
1567   InitPad();
1568   gPad->SetTopMargin(0.05);
1569   gPad->SetBottomMargin(0.15);
1570   corrX->DrawCopy();
1571   Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1572   Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
1573
1574   canvas->cd(2);
1575   InitPad();
1576   gPad->SetTopMargin(0.05);
1577   gPad->SetBottomMargin(0.15);
1578   corrY->Draw();
1579   Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1580   Tl->DrawLatex(0.5, 0.8, "|vtx-z| < 10 cm");
1581
1582   canvas->cd(3);
1583   InitPad();
1584   gPad->SetTopMargin(0.05);
1585   gPad->SetBottomMargin(0.15);
1586   gPad->SetLogx();
1587   corrZ->Draw();
1588   corrZ->GetXaxis()->SetLabelOffset(0.005);
1589   corrZ->GetXaxis()->SetTitleOffset(1.2);
1590   Tl->DrawLatex(0.5, 0.88, "|vtx-z| < 10 cm");
1591   Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
1592
1593   return canvas;
1594 }
1595
1596 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1597 {
1598   gSystem->Load("libPWG0base");
1599
1600   Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
1601
1602   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1603   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1604   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1605
1606   Prepare1DPlot(corrX);
1607   Prepare1DPlot(corrY);
1608   Prepare1DPlot(corrZ);
1609
1610   corrX->SetTitle("a) z projection");
1611   corrY->SetTitle("a) #eta projection");
1612   corrZ->SetTitle("b) p_{T} projection");
1613
1614   corrY->GetYaxis()->SetTitle("correction factor");
1615   corrZ->GetYaxis()->SetTitle("correction factor");
1616
1617   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1618
1619   TString canvasName;
1620   canvasName.Form("Track2Particle1D_%s", folder);
1621   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1622   canvas->Divide(3, 1);
1623
1624   canvas->cd(1);
1625   InitPad();
1626   corrX->DrawCopy();
1627
1628   canvas->cd(2);
1629   InitPad();
1630   corrY->Draw();
1631
1632   canvas->cd(3);
1633   InitPad();
1634   corrZ->Draw();
1635
1636   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
1637   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
1638
1639   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
1640
1641   canvasName.Form("Track2Particle1D_%s_etapt", folder);
1642   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
1643   canvas->Divide(2, 1);
1644
1645   canvas->cd(1);
1646   InitPad();
1647   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
1648   corrY->GetYaxis()->SetRangeUser(1, 1.5);
1649   corrY->GetYaxis()->SetTitleOffset(1.5);
1650   corrY->DrawCopy();
1651   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
1652   pave->AddText("|z| < 10 cm");
1653   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
1654   pave->Draw();
1655
1656   canvas->cd(2);
1657   InitPad();
1658   gPad->SetLogx();
1659   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
1660   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
1661   corrZ->GetYaxis()->SetTitleOffset(1.5);
1662   corrZ->DrawCopy();
1663   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
1664   pave->AddText("|z| < 10 cm");
1665   pave->AddText("|#eta| < 0.8");
1666   pave->Draw();
1667
1668   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
1669   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
1670 }
1671
1672 /*
1673 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1674 {
1675   gSystem->Load("libPWG0base");
1676
1677   // particle type
1678   for (Int_t particle=0; particle<4; ++particle)
1679   {
1680     TString dirName;
1681     dirName.Form("correction_%d", particle);
1682     Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
1683
1684     TString tmpx, tmpy, tmpz;
1685     tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1686     tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1687     tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
1688
1689     TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1690     TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1691     TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
1692
1693     Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
1694
1695     TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1696     TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1697     TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
1698
1699     posX->Divide(negX);
1700     posY->Divide(negY);
1701     posZ->Divide(negZ);
1702
1703     Prepare1DPlot(posX);
1704     Prepare1DPlot(posY);
1705     Prepare1DPlot(posZ);
1706
1707     Float_t min = 0.8;
1708     Float_t max = 1.2;
1709
1710     posX->SetMinimum(min);
1711     posX->SetMaximum(max);
1712     posY->SetMinimum(min);
1713     posY->SetMaximum(max);
1714     posZ->SetMinimum(min);
1715     posZ->SetMaximum(max);
1716
1717     posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1718
1719     posX->GetYaxis()->SetTitleOffset(1.7);
1720     posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1721     posY->GetYaxis()->SetTitleOffset(1.7);
1722     posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1723     posZ->GetYaxis()->SetTitleOffset(1.7);
1724     posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
1725
1726     posZ->GetXaxis()->SetRangeUser(0, 1);
1727
1728     TString canvasName;
1729     canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
1730
1731     TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1732     canvas->Divide(3, 1);
1733
1734     canvas->cd(1);
1735     InitPad();
1736     posX->DrawCopy();
1737
1738     canvas->cd(2);
1739     InitPad();
1740     posY->DrawCopy();
1741
1742     canvas->cd(3);
1743     InitPad();
1744     posZ->DrawCopy();
1745
1746     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1747     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1748   }
1749 }
1750 */
1751
1752 void CompareTrack2Particle1D(const char* file1, const char* file2, Float_t upperPtLimit = 9.9)
1753 {
1754   loadlibs();
1755
1756   const char* folderName = "dndeta_correction";
1757
1758   c = new TCanvas("CompareTrack2Particle1D", "CompareTrack2Particle1D", 1200, 400);
1759   c->Divide(3, 1);
1760
1761   for (Int_t fileId = 0; fileId < 2; fileId++)
1762   {
1763     const char* file = ((fileId == 0) ? file1 : file2);
1764     Correction1DCreatePlots(file, folderName, upperPtLimit, 1);
1765
1766     TH1* corr[3];
1767     corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1768     corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y_div_measured_y"));
1769     corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1770     /*corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x"))->Clone(Form("hist_x_%d", fileId));
1771     corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y"))->Clone(Form("hist_y_%d", fileId));
1772     corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z"))->Clone(Form("hist_z_%d", fileId));*/
1773
1774     for (Int_t i=0; i<3; i++)
1775     {
1776       c->cd(i+1);
1777       InitPad();
1778       corr[i]->GetYaxis()->SetRangeUser(0.8, 2);
1779       corr[i]->SetLineColor(fileId+1);
1780       corr[i]->DrawCopy((fileId == 0) ? "" : "SAME");
1781     }
1782   }
1783
1784   return;
1785
1786   c->SaveAs(Form("%s.gif", canvas->GetName()));
1787   c->SaveAs(Form("%s.eps", canvas->GetName()));
1788 }
1789
1790 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1791 {
1792   TFile::Open(fileName);
1793   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folder, folder);
1794   dNdEtaCorrection->LoadHistograms();
1795
1796   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1797   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1798
1799   gene->GetZaxis()->SetRangeUser(0.2, 10);
1800   meas->GetZaxis()->SetRangeUser(0.2, 10);
1801   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
1802   gene->GetZaxis()->SetRange(0, 0);
1803   meas->GetZaxis()->SetRange(0, 0);
1804
1805   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1806   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1807   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
1808   gene->GetYaxis()->SetRange(0, 0);
1809   meas->GetYaxis()->SetRange(0, 0);
1810
1811   gene->GetXaxis()->SetRangeUser(-10, 10);
1812   meas->GetXaxis()->SetRangeUser(-10, 10);
1813   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
1814   gene->GetXaxis()->SetRange(0, 0);
1815   meas->GetXaxis()->SetRange(0, 0);
1816 }
1817
1818 TCanvas* Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1819 {
1820   gSystem->Load("libPWG0base");
1821
1822   Track2Particle2DCreatePlots(fileName, folder);
1823
1824   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1825   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1826   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
1827
1828   Prepare2DPlot(corrYX);
1829   Prepare2DPlot(corrZX);
1830   Prepare2DPlot(corrZY);
1831
1832   const char* title = "";
1833   corrYX->SetTitle(title);
1834   corrZX->SetTitle(title);
1835   corrZY->SetTitle(title);
1836
1837   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1838   canvas->Divide(3, 1);
1839
1840   canvas->cd(1);
1841   InitPadCOLZ();
1842   corrYX->Draw("COLZ");
1843
1844   canvas->cd(2);
1845   InitPadCOLZ();
1846   corrZX->Draw("COLZ");
1847
1848   canvas->cd(3);
1849   InitPadCOLZ();
1850   corrZY->Draw("COLZ");
1851
1852   canvas->SaveAs(Form("corr_track2particle_%d.gif", gMax));
1853   canvas->SaveAs(Form("corr_track2particle_%d.eps", gMax));
1854
1855   return canvas;
1856 }
1857
1858 void CompareTrack2Particle2D()
1859 {
1860   gSystem->Load("libPWG0base");
1861
1862   Track2Particle2DCreatePlots("correction_maponly-positive.root");
1863
1864   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("pos_yx"));
1865   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("pos_zx"));
1866   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("pos_zy"));
1867
1868   Track2Particle2DCreatePlots("correction_maponly-negative.root");
1869
1870   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("neg_yx"));
1871   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("neg_zx"));
1872   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("neg_zy"));
1873
1874   posYX->Divide(negYX);
1875   posZX->Divide(negZX);
1876   posZY->Divide(negZY);
1877
1878   Prepare2DPlot(posYX);
1879   Prepare2DPlot(posZX);
1880   Prepare2DPlot(posZY);
1881
1882   Float_t min = 0.8;
1883   Float_t max = 1.2;
1884
1885   posYX->SetMinimum(min);
1886   posYX->SetMaximum(max);
1887   posZX->SetMinimum(min);
1888   posZX->SetMaximum(max);
1889   posZY->SetMinimum(min);
1890   posZY->SetMaximum(max);
1891
1892   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1893   canvas->Divide(3, 1);
1894
1895   canvas->cd(1);
1896   InitPadCOLZ();
1897   posYX->Draw("COLZ");
1898
1899   canvas->cd(2);
1900   InitPadCOLZ();
1901   posZX->Draw("COLZ");
1902
1903   canvas->cd(3);
1904   InitPadCOLZ();
1905   posZY->Draw("COLZ");
1906
1907   canvas->SaveAs("CompareTrack2Particle2D.gif");
1908   canvas->SaveAs("CompareTrack2Particle2D.eps");
1909 }
1910
1911 void Track2Particle3D()
1912 {
1913   // get left margin proper
1914
1915   TFile* file = TFile::Open("correction_map.root");
1916
1917   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1918
1919   corr->SetTitle("Correction Factor");
1920   SetRanges(corr->GetZaxis());
1921
1922   Prepare3DPlot(corr);
1923
1924   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1925   canvas->SetTheta(29.428);
1926   canvas->SetPhi(16.5726);
1927
1928   corr->Draw();
1929
1930   canvas->SaveAs("Track2Particle3D.gif");
1931   canvas->SaveAs("Track2Particle3D.eps");
1932 }
1933
1934 void Track2Particle3DAll()
1935 {
1936   TFile* file = TFile::Open("correction_map.root");
1937
1938   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1939   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1940   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1941
1942   gene->SetTitle("Generated Particles");
1943   meas->SetTitle("Measured Tracks");
1944   corr->SetTitle("Correction Factor");
1945
1946   Prepare3DPlot(gene);
1947   Prepare3DPlot(meas);
1948   Prepare3DPlot(corr);
1949
1950   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1951   canvas->Divide(3, 1);
1952
1953   canvas->cd(1);
1954   InitPad();
1955   gene->Draw();
1956
1957   canvas->cd(2);
1958   meas->Draw();
1959
1960   canvas->cd(3);
1961   corr->Draw();
1962
1963   canvas->SaveAs("Track2Particle3DAll.gif");
1964   canvas->SaveAs("Track2Particle3DAll.eps");
1965 }
1966
1967 void MultiplicityMC(Int_t xRangeMax = 50)
1968 {
1969   TFile* file = TFile::Open("multiplicityMC.root");
1970
1971   if (!file)
1972   {
1973     printf("multiplicityMC.root could not be opened.\n");
1974     return;
1975   }
1976
1977   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1978   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1979   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1980
1981   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1982   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1983   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1984   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1985   {
1986     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1987     proj->Fit("gaus", "0");
1988     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1989     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1990
1991     continue;
1992
1993     // draw for debugging
1994     new TCanvas;
1995     proj->DrawCopy();
1996     proj->GetFunction("gaus")->DrawCopy("SAME");
1997   }
1998
1999   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
2000
2001   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
2002   {
2003     Float_t mean = correction->GetBinContent(i);
2004     Float_t width = correctionWidth->GetBinContent(i);
2005
2006     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
2007     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
2008     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
2009
2010     for (Int_t j=fillBegin; j <= fillEnd; ++j)
2011     {
2012       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
2013     }
2014   }
2015
2016   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
2017   fMultiplicityESDCorrectedRebinned->Rebin(10);
2018   fMultiplicityESDCorrectedRebinned->Scale(0.1);
2019
2020   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
2021   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
2022   ratio->Divide(fMultiplicityMC);
2023
2024   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
2025   ratio2->Divide(fMultiplicityMC);
2026
2027   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
2028   canvas->Divide(3, 2);
2029
2030   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
2031   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
2032   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
2033   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
2034   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
2035   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
2036   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
2037
2038   canvas->cd(1); //InitPad();
2039   fMultiplicityESD->Draw();
2040   fMultiplicityMC->SetLineColor(2);
2041   fMultiplicityMC->Draw("SAME");
2042
2043   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
2044   legend->AddEntry(fMultiplicityESD, "ESD");
2045   legend->AddEntry(fMultiplicityMC, "MC");
2046   legend->Draw();
2047
2048   canvas->cd(2);
2049   fCorrelation->Draw("COLZ");
2050
2051   canvas->cd(3);
2052   correction->Draw();
2053   //correction->Fit("pol1");
2054   correctionWidth->SetLineColor(2);
2055   correctionWidth->Draw("SAME");
2056
2057   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
2058   legend->AddEntry(correction, "#bar{x}");
2059   legend->AddEntry(correctionWidth, "#sigma");
2060   legend->Draw();
2061
2062   canvas->cd(4);
2063   ratio->Draw();
2064
2065   ratio2->SetLineColor(2);
2066   ratio2->Draw("SAME");
2067
2068   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
2069   legend->AddEntry(ratio, "uncorrected");
2070   legend->AddEntry(ratio2, "corrected");
2071   legend->Draw();
2072
2073   canvas->cd(5);
2074   fMultiplicityESDCorrected->SetLineColor(kBlue);
2075   fMultiplicityESDCorrected->Draw();
2076   fMultiplicityMC->Draw("SAME");
2077   fMultiplicityESD->Draw("SAME");
2078
2079   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
2080   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
2081   legend->AddEntry(fMultiplicityMC, "MC");
2082   legend->AddEntry(fMultiplicityESD, "ESD");
2083   legend->Draw();
2084
2085   canvas->cd(6);
2086   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
2087   fMultiplicityESDCorrectedRebinned->Draw();
2088   fMultiplicityMC->Draw("SAME");
2089
2090   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
2091   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
2092   legend->AddEntry(fMultiplicityMC, "MC");
2093   legend->Draw();
2094
2095   canvas->SaveAs("MultiplicityMC.gif");
2096 }
2097
2098 void MultiplicityESD()
2099 {
2100   TFile* file = TFile::Open("multiplicityESD.root");
2101
2102   if (!file)
2103   {
2104     printf("multiplicityESD.root could not be opened.\n");
2105     return;
2106   }
2107
2108   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
2109
2110   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
2111
2112   fMultiplicityESD->Draw();
2113 }
2114
2115 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")
2116 {
2117   loadlibs();
2118
2119   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
2120   TFile::Open(correctionMapFile);
2121   dNdEtaCorrection->LoadHistograms();
2122
2123   TFile* file = TFile::Open(dataInput);
2124
2125   if (!file)
2126   {
2127     cout << "Error. File not found" << endl;
2128     return;
2129   }
2130
2131   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
2132   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
2133
2134   gROOT->cd();
2135
2136   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
2137   hist1->SetTitle("mc");
2138   Printf("mc contains %f entries", hist1->Integral());
2139   Printf("mc contains %f entries in |vtx-z| < 10, |eta| < 1, 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()));
2140
2141   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
2142   hist2->SetTitle("esd");
2143   Printf("esd contains %f entries", hist2->Integral());
2144   Printf("esd contains %f entries in |vtx-z| < 10, |eta| < 1, 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()));
2145
2146   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
2147   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
2148
2149   hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
2150   hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
2151   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
2152
2153   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
2154   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
2155   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
2156   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
2157   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
2158
2159   TH2* hist3 = (TH2*) dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram()->Clone("mc2");
2160   hist3->SetTitle("mc2");
2161   Printf("mc event contains %f entries", hist3->Integral());
2162   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()));
2163
2164   TH2* hist4 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("esd2");
2165   hist4->SetTitle("esd2");
2166   Printf("esd event contains %f entries", hist4->Integral());
2167   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()));
2168
2169   ratio = (TH2*) hist3->Clone("ratio");
2170   ratio->Divide(hist4);
2171
2172   new TCanvas; ratio->Draw("COLZ");
2173 }
2174
2175 void CompareCorrection2Generated(Float_t ptMin = 0.301, const char* dataInput = "analysis_mc.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
2176 {
2177   loadlibs();
2178
2179   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
2180   TFile::Open(correctionMapFile);
2181   dNdEtaCorrection->LoadHistograms();
2182
2183   TFile* file = TFile::Open(dataInput);
2184
2185   if (!file)
2186   {
2187     cout << "Error. File not found" << endl;
2188     return;
2189   }
2190
2191   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
2192   fdNdEtaAnalysis->LoadHistograms("dndetaTrVtx");
2193
2194   gROOT->cd();
2195
2196   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram()->Clone("mc");
2197   hist1->SetTitle("mc");
2198   Printf("mc contains %f entries", hist1->Integral());
2199   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()));
2200
2201   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram()->Clone("esd");
2202   hist2->SetTitle("esd");
2203   Printf("esd contains %f entries", hist2->Integral());
2204   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()));
2205
2206   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
2207   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
2208
2209   hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
2210   hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
2211   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
2212
2213   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
2214   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
2215   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
2216   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
2217   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
2218 }
2219
2220 void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
2221 {
2222   loadlibs();
2223
2224   TFile* file = TFile::Open(dataInput);
2225
2226   if (!file)
2227   {
2228     cout << "Error. File not found" << endl;
2229     return;
2230   }
2231
2232   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
2233   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
2234
2235   TFile* file = TFile::Open(dataInput2);
2236
2237   if (!file)
2238   {
2239     cout << "Error. File not found" << endl;
2240     return;
2241   }
2242
2243   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
2244   fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
2245
2246   gROOT->cd();
2247
2248   TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
2249   hist1->SetTitle("esd1");
2250   Printf("esd1 contains %f entries", hist1->GetEntries());
2251   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()));
2252
2253   TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
2254   hist2->SetTitle("esd2");
2255   Printf("esd2 contains %f entries", hist2->GetEntries());
2256   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()));
2257
2258   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
2259
2260   new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
2261   new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
2262   new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
2263
2264   TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
2265   TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
2266
2267   Printf("event1 contains %f entries", event1->GetEntries());
2268   Printf("event2 contains %f entries", event2->GetEntries());
2269   Printf("event1 integral is %f", event1->Integral());
2270   Printf("event2 integral is %f", event2->Integral());
2271   Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
2272   Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
2273
2274   projx1 = event1->ProjectionX();
2275   projx2 = event2->ProjectionX();
2276
2277   new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
2278
2279   projx1->Divide(projx2);
2280   new TCanvas; projx1->Draw();
2281
2282   event1->Divide(event2);
2283   new TCanvas; event1->Draw("COLZ");
2284
2285 }
2286
2287 void DrawTrackletOrigin(const char* fileName = "correction_map.root", Bool_t myFile = kTRUE)
2288 {
2289   TFile::Open(fileName);
2290
2291   Int_t maxHists = 8;
2292   TH1* hist[8];
2293
2294   const Int_t kRebin = 8;
2295
2296   const char* titles[] = { "PP", "SS", "PP'", "PS'", "PS", "SP'", "SS'", "" };
2297
2298   if (myFile)
2299   {
2300     for (Int_t i=0; i<maxHists; i++)
2301     {
2302       hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
2303       if (hist[i]->GetDimension() == 2)
2304         hist[i] = ((TH2*) hist[i])->ProjectionX(Form("fDeltaPhi_clone_%d", i));
2305     }
2306   }
2307   else
2308   {
2309     maxHists = 6;
2310     const char* names[] = { "DePhiPPTracklets", "DePhiSecTracklets", "DePhiPpTracklets", "DePhiPSTracklets", "DePhiPSdaugTracklets", "DePhiSPTracklets" };
2311     for (Int_t i=0; i<maxHists; i++)
2312       hist[i] = (TH1*) gFile->Get(names[i]);
2313   }
2314
2315   // clone before rebinning
2316   good = (TH1*) hist[0]->Clone("good");
2317   good->Add(hist[4]);
2318
2319   bad = (TH1*) hist[1]->Clone("bad");
2320   bad->Add(hist[2]);
2321   bad->Add(hist[3]);
2322   bad->Add(hist[5]);
2323   if (myFile)
2324     bad->Add(hist[6]);
2325
2326   c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c");
2327   TH1* ref = 0;
2328   Bool_t nw = kFALSE;
2329   if (!c)
2330   {
2331     c = new TCanvas("c", "c", 600, 600);
2332     nw = kTRUE;
2333     ref = (TH1*) c->GetListOfPrimitives()->At(1);
2334   }
2335   c->cd();
2336   c->SetRightMargin(0.05);
2337   c->SetTopMargin(0.05);
2338   c->SetLogy();
2339   c->SetGridx();
2340   c->SetGridy();
2341
2342   Int_t order[] = { 0, 4, 1, 2, 3, 5, 6, 7 };
2343   //Int_t colors[]  = {1,2,4,1,2,4,1,2,4};
2344   Int_t colors[]  = {1,2,3,4,6,7,8,102};
2345   Int_t markers[]  = {20, 21, 22, 23, 24, 25, 26, 27, 28};
2346
2347   TLegend* legend = new TLegend(0.75, 0.6, 0.93, 0.93);
2348   legend->SetFillColor(0);
2349   legend->SetTextSize(0.04);
2350
2351   Int_t total = 0;
2352   for (Int_t ii=0; ii<maxHists; ii++)
2353   {
2354     i = order[ii];
2355
2356     hist[i]->Rebin(kRebin);
2357     hist[i]->SetStats(kFALSE);
2358     hist[i]->SetLineColor(colors[i]);
2359     hist[i]->SetLineWidth(2);
2360     //hist[i]->SetMarkerStyle(markers[i]);
2361     //hist[i]->SetMarkerColor(colors[i]);
2362     //hist[i]->SetLineStyle(ii+1);
2363     hist[i]->GetXaxis()->SetRangeUser(-0.09, 0.09);
2364     hist[i]->GetYaxis()->SetRangeUser(5, hist[i]->GetMaximum() * 2);
2365     hist[i]->GetYaxis()->SetTitleOffset(1.3);
2366     hist[i]->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2367
2368     if (i == 0 && ref)
2369       hist[i]->Scale(1.0 / hist[i]->GetMaximum() * ref->GetMaximum());
2370
2371     hist[i]->DrawCopy(((i == 0 && nw) ? "" : "SAME"));
2372
2373     total += hist[i]->GetEntries();
2374
2375     if (i != 7)
2376       legend->AddEntry(hist[i], titles[i], "L");
2377   }
2378
2379   legend->Draw();
2380   c->SaveAs("spd_tracklets_deltaphi_detailed.eps");
2381
2382   Printf("Total: %d", total);
2383   for (Int_t i=0; i<maxHists; i++)
2384     Printf("Histogram %d (%s) contains %.2f %% of the entries", i, titles[i], 100.0 * hist[i]->GetEntries() / total);
2385
2386   printf("|  Delta phi  |  Acc. %%  |  ");
2387   for (Int_t i=0; i<maxHists; i++)
2388     printf("%3s %%   |  ", titles[i]);
2389   Printf("");
2390
2391   for (Float_t f = 0.01; f < 0.09; f += 0.01)
2392   {
2393     Int_t integralBegin = hist[0]->GetXaxis()->FindBin(-f);
2394     Int_t integralEnd = hist[0]->GetXaxis()->FindBin(f);
2395
2396     Int_t total2 = 0;
2397     for (Int_t i=0; i<maxHists; i++)
2398       total2 += (Int_t) hist[i]->Integral(integralBegin, integralEnd);
2399
2400     printf("|    %.2f     |  %6.2f  |  ", f, 100.0 * total2 / total);
2401
2402     for (Int_t i=0; i<maxHists; i++)
2403       printf("%6.2f  |  ", (hist[i]->GetEntries() > 0) ? (100.0 * hist[i]->Integral(integralBegin, integralEnd) / hist[i]->GetEntries()) : -1.0);
2404     Printf("");
2405   }
2406
2407   eff = new TH1F("eff", ";#Delta#varphi cut (rad.)", 101,-0.0005, 0.1005);
2408   cont = new TH1F("cont", "cont", 101,-0.0005, 0.1005);
2409   signalOverBg = new TH1F("signalOverBg", "signalOverBg", 101,-0.0005, 0.1005);
2410   for (Float_t cut=0.000; cut<0.10; cut += 0.001)
2411   {
2412     Float_t accGood = good->Integral(good->GetXaxis()->FindBin(-cut), good->GetXaxis()->FindBin(cut));
2413     Float_t accBad = bad->Integral(bad->GetXaxis()->FindBin(-cut), bad->GetXaxis()->FindBin(cut));
2414     Float_t sB = accGood / accBad;
2415     eff->Fill(cut, 100.0 * accGood / good->Integral());
2416     cont->Fill(cut, 100.0 * accBad / (accGood + accBad));
2417     signalOverBg->Fill(cut, sB);
2418   }
2419
2420   //new TCanvas; signalOverBg->Draw();
2421
2422   c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c2");
2423   Bool_t nw = kFALSE;
2424   if (!c)
2425   {
2426     c = new TCanvas("c2", "c2", 600, 600);
2427     nw = kTRUE;
2428   }
2429   c->cd();
2430   c->SetRightMargin(0.05);
2431   c->SetTopMargin(0.05);
2432   c->SetGridx();
2433   c->SetGridy();
2434   gPad->SetLogy();
2435   good->Rebin(kRebin);
2436   bad->Rebin(kRebin);
2437   good->GetXaxis()->SetRangeUser(-0.09, 0.09);
2438   good->GetYaxis()->SetTitleOffset(1.3);
2439   good->SetStats(0);
2440   good->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2441   good->DrawCopy((nw) ? "" : "SAME");
2442
2443   bad->SetLineColor(2);
2444   bad->SetLineStyle(2);
2445   bad->SetLineWidth(2);
2446   //bad->SetMarkerColor(2);
2447   //bad->SetMarkerStyle(7);
2448   bad->DrawCopy("SAME");
2449
2450   TLegend* legend = new TLegend(0.2, 0.13, 0.85, 0.25);
2451   legend->SetFillColor(0);
2452   legend->SetTextSize(0.04);
2453   legend->AddEntry(good, "Primaries", "L");
2454   legend->AddEntry(bad, "Secondaries + Background", "L");
2455   legend->Draw();
2456
2457   c->SaveAs("spd_tracklets_deltaphi.eps");
2458
2459   c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c3");
2460   Bool_t nw = kFALSE;
2461   if (!c)
2462   {
2463     c = new TCanvas("c3", "c3", 600, 600);
2464     nw = kTRUE;
2465   }
2466   c->cd();
2467   c->SetRightMargin(0.05);
2468   c->SetTopMargin(0.05);
2469   c->SetGridx();
2470   c->SetGridy();
2471
2472   TLegend* legend = new TLegend(0.5, 0.6, 0.93, 0.75);
2473   legend->SetFillColor(0);
2474   legend->SetTextSize(0.04);
2475   legend->AddEntry(eff, "Efficiency (%)", "L");
2476   legend->AddEntry(cont, "Contamination (%)", "L");
2477
2478   eff->SetStats(0);
2479   eff->GetXaxis()->SetRangeUser(0, 0.08);
2480   eff->GetYaxis()->SetRangeUser(1e-3, 105);
2481   eff->SetLineWidth(2);
2482   eff->DrawCopy((nw) ? "" : "SAME");
2483   cont->SetLineStyle(2);
2484   cont->SetLineWidth(2);
2485   cont->SetLineColor(2);
2486   cont->DrawCopy("SAME");
2487   legend->Draw();
2488
2489   c->SaveAs("spd_tracklets_efficiency.eps");
2490 }
2491
2492 void DrawTrackletOrigin_Compare(const char* file1, const char* file2)
2493 {
2494   DrawTrackletOrigin(file1);
2495   good1 = (TH1*) gROOT->FindObject("good")->Clone("good1");
2496   bad1 = (TH1*) gROOT->FindObject("bad")->Clone("bad1");
2497
2498   DrawTrackletOrigin(file2);
2499   good2 = (TH1*) gROOT->FindObject("good")->Clone("good2");
2500   bad2 = (TH1*) gROOT->FindObject("bad")->Clone("bad2");
2501
2502   c = new TCanvas("c4", "c4", 600, 600);
2503   c->SetRightMargin(0.05);
2504   c->SetTopMargin(0.05);
2505   c->SetGridx();
2506   c->SetGridy();
2507   gPad->SetLogy();
2508
2509   good1->Draw();
2510   bad1->SetLineColor(1);
2511   bad1->SetMarkerColor(1);
2512   bad1->Draw("SAME");
2513
2514   Float_t factor = (good1->Integral() + bad1->Integral()) / (good2->Integral() + bad2->Integral());
2515
2516   good2->Scale(factor);
2517   bad2->Scale(factor);
2518
2519   good2->SetLineColor(2);
2520   bad2->SetMarkerColor(2);
2521
2522   good2->Draw("SAME");
2523   bad2->Draw("SAME");
2524
2525   good1->GetYaxis()->SetRangeUser(1, TMath::Max(good1->GetMaximum(), good2->GetMaximum()) * 1.1);
2526 }
2527
2528 void Tracklets_Asymmetry()
2529 {
2530   TFile::Open("correction_map.root");
2531
2532   Int_t maxHists = 7;
2533   TH1* hist[8];
2534
2535   Int_t colors[]  = {1,2,3,4,6,7,8,102};
2536   const char* titles[] = { "PP", "SS", "PP'", "PS'", "PS", "SP'", "SS'", "" };
2537
2538   TLegend* legend = new TLegend(0.75, 0.6, 0.93, 0.93);
2539
2540   for (Int_t i=0; i<maxHists; i++)
2541   {
2542     hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
2543     hist[i]->Rebin(10);
2544
2545     for (Int_t j=hist[i]->GetNbinsX()/2; j<=hist[i]->GetNbinsX(); j++)
2546       if (hist[i]->GetBinContent(j) > 0)
2547         hist[i]->SetBinContent(j, (hist[i]->GetBinContent(j) -  hist[i]->GetBinContent(hist[i]->GetXaxis()->FindBin(-hist[i]->GetXaxis()->GetBinCenter(j)))) / hist[i]->GetBinContent(j));
2548
2549     hist[i]->SetStats(kFALSE);
2550     hist[i]->SetLineColor(colors[i]);
2551     hist[i]->GetXaxis()->SetRangeUser(0.001, 0.09);
2552     //hist[i]->GetYaxis()->SetRangeUser(5, hist[i]->GetMaximum() * 2);
2553     hist[i]->GetYaxis()->SetTitleOffset(1.3);
2554     hist[i]->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2555     hist[i]->Draw(((i == 0) ? "" : "SAME"));
2556
2557     legend->AddEntry(hist[i], titles[i], "L");
2558   }
2559
2560   legend->Draw();
2561 }
2562
2563 TH2* GetCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2564 {
2565   // returns the correction factor with pt integrated out
2566
2567   loadlibs();
2568
2569   TFile::Open(fileName);
2570
2571   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
2572   if (!dNdEtaCorrection->LoadHistograms())
2573     return;
2574
2575   //  hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
2576
2577   gener = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
2578   measu = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
2579
2580   gener->GetZaxis()->SetRange(gener->GetZaxis()->FindBin(ptmin), gener->GetNbinsZ()+1);
2581   TH2D *gener_xy = gener->Project3D("yx");
2582
2583   measu->GetZaxis()->SetRange(measu->GetZaxis()->FindBin(ptmin), measu->GetNbinsZ()+1);
2584   TH2D *measu_xy = measu->Project3D("yx");
2585
2586   cout << measu->GetZaxis()->FindBin(ptmin) << " " << measu->GetNbinsZ()+1 << endl;
2587
2588   TCanvas *canp = new TCanvas("canp","canp",600,1000);
2589   canp->Divide(1,2,0.0001,0.0001);
2590   canp->cd(1);
2591   gener_xy->Draw("COLZ");
2592   canp->cd(2);
2593   measu_xy->Draw("COLZ");
2594
2595
2596   TCanvas *canpr = new TCanvas("canpr","canpr",700,500);
2597   canpr->cd();
2598   TH2D *proj = new TH2D(*gener_xy);
2599   proj->Divide(measu_xy);
2600
2601 //   proj = hist->Project3D("yx");
2602   proj->Draw("COLZ");
2603
2604   return proj;
2605 }
2606
2607 void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2608 {
2609   TH2* proj = GetCorrection(fileName, dirName, ptmin);
2610
2611   const Float_t limit = 5;
2612
2613   TString array = "{";
2614   TString arrayEnd = "}";
2615
2616   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2617   {
2618     Int_t begin = -1;
2619     Int_t end = -1;
2620     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2621     {
2622       if (begin == -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2623         begin = x;
2624       if (begin != -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2625         end = x;
2626     }
2627     Printf("Limits for y = %d are %d to %d", y, begin, end);
2628
2629     if (y > 1)
2630       array += ", ";
2631     array += Form("%d", begin);
2632
2633     if (y > 1)
2634       arrayEnd.Prepend(", ");
2635     arrayEnd.Prepend(Form("%d", (end == -1) ? -1 : proj->GetNbinsX() + 1 - end));
2636   }
2637   array += "}";
2638   arrayEnd.Prepend("{");
2639
2640   Printf("Begin array:");
2641   Printf("%s", array.Data());
2642
2643   Printf("End array (mirrored) (should be the same):");
2644   Printf("%s", arrayEnd.Data());
2645 }
2646
2647 void AverageMultiplicity(const char* fileName = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
2648 {
2649   loadlibs();
2650
2651   TFile::Open(fileName);
2652
2653   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
2654   dNdEtaCorrection->LoadHistograms();
2655   TH2* events = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram();
2656   TH3* tracks = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram();
2657
2658   Float_t nEvents = events->Integral(events->GetXaxis()->FindBin(-1), events->GetXaxis()->FindBin(1), 0, events->GetNbinsY()+1);
2659   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);
2660
2661   Printf("%f %f --> %f", nEvents, nTracks, nTracks / nEvents);
2662 }
2663
2664 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")
2665 {
2666   loadlibs();
2667
2668   TFile::Open(rawFile);
2669   dNdEtaAnalysis* raw = new dNdEtaAnalysis("dndeta", "dndeta");
2670   raw->LoadHistograms("fdNdEtaAnalysisESD");
2671   raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->GetXaxis()->SetRangeUser(-vertexRange, vertexRange);
2672   tracks = raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Project3D("y");
2673   events = raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("events", 0, raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->GetNbinsY() + 1);
2674   Float_t nEvents = events->Integral(events->FindBin(-vertexRange), events->FindBin(vertexRange));
2675   tracks->Scale(1.0 / nEvents / tracks->GetBinWidth(1));
2676
2677   TFile::Open(mcFile);
2678   dNdEtaAnalysis* mc = new dNdEtaAnalysis("dndeta", "dndeta");
2679   mc->LoadHistograms("dndetaTrVtx");
2680   mcH = mc->GetdNdEtaPtCutOffCorrectedHistogram(0);
2681
2682   new TCanvas;
2683   mcH->SetLineColor(2);
2684   mcH->DrawCopy();
2685   tracks->DrawCopy("SAME");
2686
2687   new TCanvas;
2688   mcH->GetYaxis()->SetRangeUser(0, 5);
2689   mcH->Divide(tracks);
2690   mcH->DrawCopy();
2691   mcH->Fit("pol0", "", "", -etaRange, etaRange);
2692 }
2693
2694 void TrackCuts_Comparison_MC(char* histName, Int_t plotWhich = 0, const char* fileName = "correction_map.root", Bool_t mirror = kFALSE)
2695 {
2696   // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
2697   //    --> manually disable it in the run.C
2698   //
2699   // plotWhich: 0 = only before
2700   //            1 = both
2701   //            2 = only after
2702   //
2703   // mirror: kTRUE --> project negative values on the positive side
2704
2705
2706   file = TFile::Open(fileName);
2707
2708   Int_t count = 0;
2709   Int_t colors[] = { 1, 2, 3, 4, 5, 6 };
2710
2711   TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
2712   TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
2713   TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
2714
2715   TCanvas* c1 = new TCanvas(histName, histName, 800, 1200);
2716   c1->Divide(1, 2);
2717   //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
2718   //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
2719
2720   const char* folders2[] = { "before_cuts", "after_cuts" };
2721   Bool_t first = kTRUE;
2722   for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
2723   {
2724     const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
2725     const char* names[] =    { "all", "primaries", "secondaries" };
2726     TH1* base = 0;
2727     TH1* prim = 0;
2728     TH1* sec = 0;
2729     for (Int_t i = 0; i < 3; i++)
2730     {
2731       TString folder;
2732       folder.Form("%s/%s/%s", folders1[i], folders2[j], histName);
2733       TH1* hist = (TH1*) file->Get(folder);
2734
2735       if (mirror)
2736       {
2737         for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
2738         {
2739           Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
2740           if (bin != newBin)
2741           {
2742             hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
2743             hist->SetBinContent(bin, 0);
2744           }
2745         }
2746       }
2747
2748       legend->AddEntry(hist, Form("%s %s", names[i], folders2[j]));
2749
2750       c1->cd(1);
2751       hist->SetLineColor(colors[count]);
2752       hist->DrawCopy((count == 0) ? "" : "SAME");
2753
2754       switch (i)
2755       {
2756         case 0: base = hist; break;
2757         case 1: prim = hist; break;
2758         case 2: sec = hist; break;
2759       }
2760
2761       count++;
2762     }
2763
2764     TH1* eff    = (TH1*) prim->Clone("eff"); eff->Reset();
2765     TH1* purity = (TH1*) prim->Clone("purity"); purity->Reset();
2766
2767     for (Int_t bin = 1; bin <= prim->GetNbinsX(); bin++)
2768     {
2769       eff->SetBinContent(bin, prim->Integral(1, bin) / prim->Integral(1, prim->GetNbinsX() + 1));
2770       if (prim->Integral(1, bin) + sec->Integral(1, bin) > 0)
2771         purity->SetBinContent(bin, sec->Integral(1, bin) / (prim->Integral(1, bin) + sec->Integral(1, bin)));
2772     }
2773
2774     eff->GetYaxis()->SetRangeUser(0, 1);
2775     eff->SetLineColor(colors[0+j*2]);
2776     eff->SetStats(kFALSE);
2777     purity->SetLineColor(colors[1+j*2]);
2778
2779     legend3->AddEntry(eff, Form("%s: efficiency", folders2[j]));
2780     legend3->AddEntry(purity, Form("%s: contamination", folders2[j]));
2781
2782     c1->cd(2);
2783     eff->DrawCopy((first) ? "" : "SAME");
2784     first = kFALSE;
2785     purity->DrawCopy("SAME");
2786   }
2787
2788   c1->cd(1)->SetLogy();
2789   c1->cd(1)->SetGridx();
2790   c1->cd(1)->SetGridy();
2791   legend->Draw();
2792
2793   //c2->cd();
2794  // c2->SetGridx();
2795  // c2->SetGridy();
2796   //legend2->Draw();
2797
2798   c1->cd(2)->SetGridx();
2799   c1->cd(2)->SetGridy();
2800   legend3->Draw();
2801
2802   //c1->SaveAs(Form("%s.png", histName));
2803 }
2804
2805 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")
2806 {
2807   // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
2808   //    --> manually disable it in the run.C
2809   //
2810   // plotWhich: 0 = only before
2811   //            1 = both
2812   //            2 = only after
2813   //
2814   // mirror: kTRUE --> project negative values on the positive side
2815
2816
2817   Int_t count = 0;
2818   Int_t colors[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
2819
2820   TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
2821   legend->SetFillColor(0);
2822   TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
2823   TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
2824
2825   TCanvas* c1 = new TCanvas(histName, histName, 600, 600);
2826   //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
2827   //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
2828
2829   const char* folders2[] = { "before_cuts", "after_cuts" };
2830   Bool_t first = kTRUE;
2831   for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
2832   {
2833     const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
2834     const char* names[] =    { "all", "primaries", "secondaries" };
2835
2836     Float_t normalize[3];
2837
2838     for (Int_t i = 0; i < 2; i++)
2839     {
2840       file = TFile::Open((i == 0) ? fileName1 : fileName2);
2841
2842       for (Int_t k = 1; k < 3; k++)
2843       {
2844         TString folder;
2845         folder.Form("%s/%s/%s", folders1[k], folders2[j], histName);
2846         Printf("%s", folder.Data());
2847         TH1* hist = (TH1*) file->Get(folder);
2848
2849         if (mirror)
2850         {
2851           for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
2852           {
2853             Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
2854             if (bin != newBin)
2855             {
2856               hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
2857               hist->SetBinContent(bin, 0);
2858             }
2859           }
2860         }
2861
2862         if (i == 0)
2863         {
2864           normalize[k] = hist->Integral();
2865         }
2866         else
2867           hist->Scale(normalize[k] / hist->Integral());
2868
2869         legend->AddEntry(hist, Form("%s %s %s", (i == 0) ? label1 : label2, (k == 1) ? "primaries" : "secondaries", folders2[j]));
2870
2871         c1->cd();
2872         hist->SetStats(0);
2873         hist->SetLineColor(colors[count]);
2874         hist->DrawCopy((count == 0) ? "" : "SAME");
2875
2876         count++;
2877       }
2878     }
2879
2880   }
2881
2882   //c1->SetLogy();
2883   c1->SetGridx();
2884   c1->SetGridy();
2885   legend->Draw();
2886 }
2887
2888 void TrackCuts_DCA()
2889 {
2890   file = TFile::Open("correction_map.root");
2891   hist = (TH2*) file->Get("esd_track_cuts/before_cuts/dXYvsDZ");
2892
2893   TCanvas* c1 = new TCanvas("c1", "c1", 600, 600);
2894   c1->SetLogz();
2895   c1->SetRightMargin(0.12);
2896   c1->SetBottomMargin(0.12);
2897
2898   hist->SetStats(kFALSE);
2899   hist->Draw("COLZ");
2900
2901   ellipse = new TEllipse(0, 0, 4);
2902   ellipse->SetLineWidth(2);
2903   ellipse->SetLineStyle(2);
2904   ellipse->SetFillStyle(0);
2905   ellipse->Draw();
2906
2907   c1->SaveAs("trackcuts_dca_2d.eps");
2908 }
2909
2910 void FindNSigma(TH2* hist, Int_t nSigma = 1)
2911 {
2912   TH1* proj = hist->ProjectionY();
2913   proj->Reset();
2914
2915   for (Int_t bin=1; bin<=proj->GetNbinsX(); bin++)
2916   {
2917     if (hist->Integral(1, hist->GetNbinsX(), bin, bin) == 0)
2918       continue;
2919
2920     Int_t limit = -1;
2921     for (limit = 1; limit<=hist->GetNbinsX(); limit++)
2922   }
2923 }
2924
2925 void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2926 {
2927   TH2* proj = GetCorrection(fileName, dirName, ptmin);
2928
2929   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2930     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2931       if (proj->GetBinContent(x, y) > 5 || proj->GetBinContent(x, y) == 0)
2932       {
2933         proj->SetBinContent(x, y, 0);
2934       }
2935       else
2936         proj->SetBinContent(x, y, 1);
2937
2938
2939   input->Multiply(proj);
2940 }
2941
2942 void MakeGaussianProfile(const char* histName = "fVertexCorrelation", Bool_t subtractMean = kFALSE)
2943 {
2944     TFile::Open("correction_map.root");
2945
2946     TH2* hist2d = (TH2*) gFile->Get(histName);
2947     hist2d->Sumw2();
2948
2949     TH1* result = hist2d->ProjectionX("result");
2950     result->GetYaxis()->SetTitle(hist2d->GetYaxis()->GetTitle());
2951     result->Reset();
2952
2953     for (Int_t x=1; x<hist2d->GetNbinsX(); ++x)
2954     {
2955         hist = hist2d->ProjectionY(Form("temp_%d", x), x, x);
2956         if (hist->GetEntries() == 0)
2957             continue;
2958         if (hist->Fit("gaus") == 0)
2959         {
2960             func = hist->GetFunction("gaus");
2961             mean = func->GetParameter(1);
2962             error = func->GetParError(1);
2963
2964             if (subtractMean)
2965                 mean = hist2d->GetXaxis()->GetBinCenter(x) - mean;
2966
2967             result->SetBinContent(x, mean);
2968             result->SetBinError(x, error);
2969
2970             if (x % 10 == 0)
2971             {
2972                 new TCanvas;
2973                 ((TH1*) hist->Clone())->DrawCopy();
2974             }
2975         }
2976         //break;
2977     }
2978
2979     new TCanvas;
2980     result->GetYaxis()->SetRangeUser(-0.2, 0.2);
2981     result->Draw();
2982 }
2983
2984 TH2* GetAcceptance(void* corr2d_void)
2985 {
2986         corr2d = (AliCorrectionMatrix2D*) corr2d_void;
2987         corr_xy = (TH2*) corr2d->GetCorrectionHistogram()->Clone("acceptance");
2988
2989         // fold in acceptance
2990         for (Int_t x=1; x<=corr_xy->GetNbinsX(); ++x)
2991                 for (Int_t y=1; y<=corr_xy->GetNbinsY(); ++y)
2992                 {
2993                         if (corr_xy->GetBinContent(x, y) > 1.5)
2994                                 corr_xy->SetBinContent(x, y, 0);
2995
2996                         if (corr_xy->GetBinContent(x, y) > 0)
2997                                 corr_xy->SetBinContent(x, y, 1);
2998
2999                         corr_xy->SetBinError(x, y, 0);
3000                 }
3001
3002         return corr_xy;
3003 }
3004
3005 void ZeroOutsideAcceptance(TH2* acc, TH3* hist)
3006 {
3007   for (Int_t x=0; x<=acc->GetNbinsX()+1; ++x)
3008     for (Int_t y=0; y<=acc->GetNbinsY()+1; ++y)
3009     {
3010       if (acc->GetBinContent(x, y) > 2 || acc->GetBinContent(x, y) == 0)
3011       {
3012         for (Int_t z=0; z<=hist->GetNbinsZ()+1; ++z)
3013         {
3014           hist->SetBinContent(x, y, z, 0);
3015           hist->SetBinError(x, y, z, 0);
3016         }
3017       }
3018     }
3019 }
3020
3021 void DrawPhi()
3022 {
3023   loadlibs();
3024
3025   TFile::Open("correction_map.root");
3026   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
3027   if (!dNdEtaCorrection->LoadHistograms())
3028     return 0;
3029
3030   TFile::Open("analysis_esd.root");
3031   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
3032   fdNdEtaAnalysis->LoadHistograms();
3033
3034   // acc. map!
3035   //acc = GetAcceptance(dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000));
3036   acc = dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000)->GetCorrectionHistogram();
3037   //new TCanvas; acc->Draw("COLZ");
3038
3039   histG = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
3040   ZeroOutsideAcceptance(acc, histG);
3041   //new TCanvas; histG->Project3D("yx")->Draw("COLZ");
3042   //histG->GetYaxis()->SetRangeUser(-0.9, 0.9);
3043
3044   histM = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
3045   ZeroOutsideAcceptance(acc, histM);
3046   //histM->GetYaxis()->SetRangeUser(-0.9, 0.9);
3047
3048   TFile::Open("analysis_mc.root");
3049   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndetaTrVtxMC", "dndetaTrVtxMC");
3050   fdNdEtaAnalysis2->LoadHistograms("dndetaTrVtx");
3051
3052   histMC = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
3053   ZeroOutsideAcceptance(acc, histMC);
3054   //new TCanvas; histMC->Project3D("yx2")->Draw("COLZ");
3055
3056   //histG->GetZaxis()->SetRangeUser(1,2); histMC->GetZaxis()->SetRangeUser(1,2);
3057   new TCanvas; a = histG->Project3D("yx3"); a->Add(histMC->Project3D("yx4"), -1); a->Draw("COLZ");
3058
3059   //histMC->GetYaxis()->SetRangeUser(-0.9, 0.9);
3060
3061   c = new TCanvas;
3062
3063   histG->GetXaxis()->SetRangeUser(-9.9, 9.9);
3064   histG->Project3D("z")->DrawCopy();
3065
3066   histM->GetXaxis()->SetRangeUser(-9.9, 9.9);
3067   proj = histM->Project3D("z2");
3068   proj->SetLineColor(2);
3069   proj->DrawCopy("SAME");
3070
3071   histMC->GetXaxis()->SetRangeUser(-9.9, 9.9);
3072   projMC = histMC->Project3D("z3");
3073   projMC->SetLineColor(4);
3074   projMC->DrawCopy("SAME");
3075 }
3076
3077 void PrintEventStats(Int_t corrID = 3, const char* fileName = "correction_map.root", const char* dir = "dndeta_correction")
3078 {
3079   loadlibs();
3080
3081   /*
3082   TFile::Open("analysis_mc.root");
3083   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD");
3084   fdNdEtaAnalysis->LoadHistograms();
3085   trackHist = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
3086   eventHist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetGeneratedHistogram();
3087   */
3088
3089   TFile::Open(fileName);
3090   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dir, dir);
3091   if (!dNdEtaCorrection->LoadHistograms())
3092     return;
3093   trackHist = dNdEtaCorrection->GetCorrection(corrID)->GetTrackCorrection()->GetGeneratedHistogram();
3094   eventHist = dNdEtaCorrection->GetCorrection(corrID)->GetEventCorrection()->GetGeneratedHistogram();
3095
3096   trackHist->GetXaxis()->SetRange(0, trackHist->GetNbinsX()+1);
3097   trackHist->GetZaxis()->SetRange(0, trackHist->GetNbinsZ()+1);
3098   eta = trackHist->Project3D("y");
3099
3100   events = eventHist->Integral(0, eventHist->GetNbinsX()+1, 0, eventHist->GetNbinsY()+1);
3101
3102   eta->Scale(1.0 / events);
3103
3104   Float_t avgN = eta->Integral(0, eta->GetNbinsX()+1);
3105   Printf("<N> = %f", avgN);
3106
3107   eta->Scale(1.0 / eta->GetXaxis()->GetBinWidth(1));
3108
3109   Printf("dndeta | eta = 0 is %f", (eta->GetBinContent(eta->FindBin(0.01)) + eta->GetBinContent(eta->FindBin(-0.01))) / 2);
3110   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));
3111   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));
3112   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));
3113
3114   stats = (TH2*) gFile->Get("fEventStats");
3115   proj = stats->ProjectionX();
3116   gROOT->ProcessLine(".L PrintHist.C");
3117   PrintHist2D(stats);
3118   PrintHist(proj);
3119
3120   Float_t ua5_SD = 0.153;
3121   Float_t ua5_DD = 0.080;
3122   Float_t ua5_ND = 0.767;
3123
3124   Printf("+++ FRACTIONS +++");
3125
3126   Printf("ND: %f", proj->GetBinContent(3) / proj->GetBinContent(1));
3127   Printf("SD: %f", proj->GetBinContent(4) / proj->GetBinContent(1));
3128   Printf("DD: %f", proj->GetBinContent(5) / proj->GetBinContent(1));
3129
3130   Printf("+++ TRIGGER EFFICIENCIES +++");
3131
3132   Printf("INEL = %.1f", 100. * (proj->GetBinContent(1) - stats->GetBinContent(1, 1) - stats->GetBinContent(1, 3)) / proj->GetBinContent(1));
3133   Printf("NSD  = %.1f", 100. * (proj->GetBinContent(2) - stats->GetBinContent(2, 1) - stats->GetBinContent(2, 3)) / proj->GetBinContent(2));
3134
3135
3136   Float_t trigND = 100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1) - stats->GetBinContent(3, 3)) / proj->GetBinContent(3);
3137   Float_t trigSD = 100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1) - stats->GetBinContent(4, 3)) / proj->GetBinContent(4);
3138   Float_t trigDD = 100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1) - stats->GetBinContent(5, 3)) / proj->GetBinContent(5);
3139
3140   Printf("ND  = %.1f", trigND);
3141   Printf("SD  = %.1f", trigSD);
3142   Printf("DD  = %.1f", trigDD);
3143
3144   Float_t trigINELUA5 = ua5_SD * trigSD + ua5_DD * trigDD + ua5_ND * trigND;
3145   Float_t trigNSDUA5  = (ua5_DD * trigDD + ua5_ND * trigND) / (ua5_DD + ua5_ND);
3146   Printf("INEL (UA5)  = %.1f", trigINELUA5);
3147   Printf("NSD (UA5)  = %.1f", trigNSDUA5);
3148
3149   Printf("+++ VERTEX EFFICIENCIES +++");
3150
3151   Printf("INEL = %.1f", 100. * (stats->GetBinContent(1, 3) + stats->GetBinContent(1, 4)) / proj->GetBinContent(1));
3152   Printf("NSD  = %.1f", 100. * (stats->GetBinContent(2, 3) + stats->GetBinContent(2, 4)) / proj->GetBinContent(2));
3153
3154   Float_t vtxND = 100. * (stats->GetBinContent(3, 3) + stats->GetBinContent(3, 4)) / proj->GetBinContent(3);
3155   Float_t vtxSD = 100. * (stats->GetBinContent(4, 3) + stats->GetBinContent(4, 4)) / proj->GetBinContent(4);
3156   Float_t vtxDD = 100. * (stats->GetBinContent(5, 3) + stats->GetBinContent(5, 4)) / proj->GetBinContent(5);
3157   Printf("ND  = %.1f", vtxND);
3158   Printf("SD  = %.1f", vtxSD);
3159   Printf("DD  = %.1f", vtxDD);
3160
3161   Float_t vtxINELUA5 = ua5_SD * vtxSD + ua5_DD * vtxDD + ua5_ND * vtxND;
3162   Float_t vtxNSDUA5  = (ua5_DD * vtxDD + ua5_ND * vtxND) / (ua5_DD + ua5_ND);
3163   Printf("INEL (UA5)  = %.1f", vtxINELUA5);
3164   Printf("NSD (UA5)  = %.1f", vtxNSDUA5);
3165
3166   Printf("+++ TRIGGER + VERTEX EFFICIENCIES +++");
3167
3168   Printf("INEL = %.1f", 100. * stats->GetBinContent(1, 4) / proj->GetBinContent(1));
3169   Printf("NSD  = %.1f", 100. * stats->GetBinContent(2, 4) / proj->GetBinContent(2));
3170   Printf("ND  = %.1f",  100. * stats->GetBinContent(3, 4) / proj->GetBinContent(3));
3171   Printf("SD  = %.1f",  100. * stats->GetBinContent(4, 4) / proj->GetBinContent(4));
3172   Printf("DD  = %.1f",  100. * stats->GetBinContent(5, 4) / proj->GetBinContent(5));
3173
3174
3175
3176   for (Int_t i=7; i<=proj->GetNbinsX(); i++)
3177     if (proj->GetBinContent(i) > 0)
3178       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));
3179
3180   //eta->Draw();
3181 }
3182
3183 void TestAsymmetry()
3184 {
3185   loadlibs();
3186
3187   TFile* file2 = TFile::Open("analysis_mc.root");
3188
3189   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
3190   fdNdEtaAnalysis->LoadHistograms();
3191   fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
3192
3193   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy();
3194
3195   hist = (TH1*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
3196   hist2 = (TH1*) hist->Clone("hist2");
3197   for (Int_t x=1; x<=hist->GetNbinsX(); x++)
3198     for (Int_t y=1; y<=hist->GetNbinsY(); y++)
3199       for (Int_t z=1; z<=hist->GetNbinsZ(); z++)
3200       {
3201         Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
3202         hist->SetBinContent(x, y, z, hist2->GetBinContent(hist->GetNbinsX() / 2, TMath::Min(y, hist->GetNbinsY() + 1 - y), z));
3203       }
3204
3205   hist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
3206   for (Int_t x=1; x<=hist->GetNbinsX(); x++)
3207     for (Int_t y=1; y<=hist->GetNbinsY(); y++)
3208       {
3209         //Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
3210         hist->SetBinContent(x, y, hist->GetBinContent(hist->GetNbinsX() / 2, y));
3211       }
3212
3213   fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
3214   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerColor(2);
3215   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetLineColor(2);
3216   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerStyle(5);
3217   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy("SAMEP");
3218 }
3219
3220 void DeltaPhiFromPt(Float_t smearing = 0.005)
3221 {
3222   loadlibs();
3223
3224   TFile::Open("analysis_mc.root");
3225   hist = (TH1*) gFile->Get("dndeta_check_pt");
3226
3227   dPhiHist = new TH1F("dPhiHist", ";#Delta phi", 400, -0.1, 0.1);
3228   dPhiHist2 = new TH1F("dPhiHist2", ";#Delta phi", 400, -0.1, 0.1);
3229
3230   for (Int_t i=1; i<=hist->GetNbinsX(); i++)
3231   {
3232     Float_t pt = hist->GetBinCenter(i);
3233     Float_t deltaPhi = (0.076 - 0.039) / 2 / (pt / 0.15);
3234
3235     if (smearing > 0)
3236     {
3237       gaus = new TF1("mygaus", "gaus(0)", -0.1, 0.1);
3238       gaus->SetParameters(1, -deltaPhi, smearing);
3239
3240       dPhiHist->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
3241
3242       dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
3243       gaus->SetParameters(1, deltaPhi, smearing);
3244       dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
3245     }
3246     else
3247 {
3248 dPhiHist->Fill(deltaPhi, hist->GetBinContent(i) / 2);
3249 dPhiHist2->Fill(deltaPhi, hist->GetBinContent(i) / 2);
3250 dPhiHist2->Fill(-deltaPhi, hist->GetBinContent(i) / 2);
3251 }
3252   }
3253
3254   new TCanvas;
3255   dPhiHist->Draw();
3256   dPhiHist2->SetLineColor(2);
3257   dPhiHist2->Draw("SAME");
3258   gPad->SetLogy();
3259
3260   TFile::Open("trackletsDePhi.root");
3261   //TFile::Open("tmp/correction_maponly-positive.root");
3262   //TFile::Open("tmp/correction_map.root");
3263   //tracklets = (TH1*) gFile->Get(Form("fDeltaPhi_%d", 0));
3264   tracklets = (TH1*) gFile->Get("DePhiPPTracklets");
3265   tracklets->Scale(1.0 / tracklets->GetMaximum() * dPhiHist->GetMaximum());
3266   tracklets->SetLineColor(4);
3267   tracklets->Draw("SAME");
3268 }
3269
3270 void VertexDistributions()
3271 {
3272   loadlibs();
3273
3274   TFile::Open("correction_map.root");
3275   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
3276   if (!dNdEtaCorrection->LoadHistograms())
3277     return;
3278
3279   all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all");
3280   trigger = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("trigger");
3281   vtx = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("vtx");
3282
3283   nottriggered = (TH1*) all->Clone("nottriggered");
3284   nottriggered->Add(trigger, -1);
3285
3286   novertex = (TH1*) trigger->Clone("novertex");
3287   novertex->Add(vtx, -1);
3288
3289   temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
3290   highmult = temphist->ProjectionX("highmult", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
3291   //all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
3292
3293   for (Int_t i=1; i<=trigger->GetNbinsX(); i++)
3294   {
3295     all->SetBinContent(i, all->GetBinContent(i) / all->GetBinWidth(i));
3296     trigger->SetBinContent(i, trigger->GetBinContent(i) / trigger->GetBinWidth(i));
3297     vtx->SetBinContent(i, vtx->GetBinContent(i) / vtx->GetBinWidth(i));
3298     nottriggered->SetBinContent(i, nottriggered->GetBinContent(i) / nottriggered->GetBinWidth(i));
3299     novertex->SetBinContent(i, novertex->GetBinContent(i) / novertex->GetBinWidth(i));
3300     highmult->SetBinContent(i, highmult->GetBinContent(i) / highmult->GetBinWidth(i));
3301   }
3302
3303   new TCanvas;
3304   vtx->SetTitle("");
3305   vtx->SetStats(0);
3306   vtx->DrawCopy("HIST");
3307
3308   all->Scale(1.0 / all->Integral());
3309   nottriggered->Scale(1.0 / nottriggered->Integral());
3310   novertex->Scale(1.0 / novertex->Integral());
3311   highmult->Scale(1.0 / highmult->Integral());
3312
3313   new TCanvas;
3314   all->Draw("HIST");
3315   novertex->SetLineColor(2);
3316   novertex->Draw("HISTSAME");
3317   highmult->SetLineColor(3);
3318   highmult->Draw("HISTSAME");
3319
3320   legend = new TLegend(0.5, 0.5, 0.8, 0.8);
3321   legend->SetFillColor(0);
3322   legend->AddEntry(all, "all");
3323   legend->AddEntry(novertex, "no vertex");
3324   legend->AddEntry(highmult, "mult > 10");
3325   legend->Draw();
3326
3327   new TCanvas;
3328   trigger->Scale(1.0 / trigger->Integral());
3329   vtx->Scale(1.0 / vtx->Integral());
3330
3331   trigger->Divide(vtx);
3332
3333   trigger->Draw();
3334   //vtx->SetLineColor(2);
3335   //vtx->Draw("SAME");
3336
3337   //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
3338   temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram();
3339   //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram();
3340
3341   temphist = (TH2*) gFile->Get("fTemp1");
3342
3343   new TCanvas;
3344   legend = new TLegend(0.7, 0.7, 0.99, 0.99);
3345   legend->SetFillColor(0);
3346
3347   Bool_t first = kTRUE;
3348   for (Int_t i=0; i<20; i+=5)
3349   {
3350     highmult = temphist->ProjectionX("highmult", i+1, i+1+4);
3351     highmult->Rebin(10);
3352     Printf("%f", highmult->Integral());
3353     if (highmult->Integral() <= 0)
3354       continue;
3355
3356     for (Int_t j=1; j<=trigger->GetNbinsX(); j++)
3357       highmult->SetBinContent(j, highmult->GetBinContent(j) / highmult->GetBinWidth(j));
3358
3359     highmult->Scale(1.0 / highmult->Integral());
3360     highmult->SetLineColor((i/5)+1);
3361     highmult->GetYaxis()->SetRangeUser(0, 0.15);
3362     if (first)
3363     {
3364       highmult->DrawCopy();
3365       first = kFALSE;
3366     }
3367     else
3368       highmult->DrawCopy("SAME");
3369     legend->AddEntry(highmult->Clone(), Form("%d <= N <= %d", i, i+4));
3370   }
3371   legend->Draw();
3372
3373 }
3374
3375 void PlotPt1DCorrection()
3376 {
3377   const char* files[] = { "field.root", "field_onlyprim.root", "nofield.root", "nofield_onlyprim.root" };
3378   const char* names[] = { "B: all", "B: primaries", "No B: all", "No B: primaries" };
3379   Int_t colors[] = { 1, 2, 3, 4 };
3380
3381   loadlibs();
3382
3383   dummy = new TH2F("dummy", ";p_{T};correction", 100, 0, 1.4, 100, 0.5, 3);
3384   dummy->SetStats(0);
3385   //dummy->GetYaxis()->SetTitleOffset(1.3);
3386   dummy->Draw();
3387
3388   legend = new TLegend(0.48, 0.57, 0.88, 0.88);
3389   legend->SetFillColor(0);
3390
3391   for (Int_t i=0; i<4; i++)
3392   {
3393     TFile::Open(files[i]);
3394     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
3395     if (!dNdEtaCorrection->LoadHistograms())
3396       return;
3397
3398     hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get1DCorrectionHistogram("z", -9.9, 9.9, -0.79, 0.79);
3399     hist->SetLineColor(colors[i]);
3400     hist->SetLineWidth(2);
3401     hist->SetMarkerColor(colors[i]);
3402     hist->Draw("SAME");
3403
3404     legend->AddEntry(hist, names[i], "L");
3405   }
3406
3407   legend->Draw();
3408 }
3409
3410 void FitDiamond()
3411 {
3412   TFile::Open("analysis_esd_raw.root");
3413
3414   hist = (TH3*) gFile->Get("vertex_check");
3415
3416   gStyle->SetOptFit(1);
3417
3418   TH1* proj[3];
3419   proj[0] = hist->ProjectionX();
3420   proj[1] = hist->ProjectionY();
3421   proj[2] = hist->ProjectionZ();
3422
3423   for (Int_t i=0; i<3; i++)
3424   {
3425     c = new TCanvas;
3426     proj[i]->Draw();
3427     proj[i]->Fit("gaus");
3428
3429     c->SaveAs(Form("FitDiamond_%d.png", i));
3430   }
3431 }
3432
3433 void CompareDiamond(const char* mc, const char* data)
3434 {
3435   TFile::Open(mc);
3436
3437   hist = (TH3*) gFile->Get("vertex_check");
3438
3439   gStyle->SetOptFit(1);
3440
3441   TH1* proj[3];
3442   proj[0] = hist->ProjectionX("vertex_check_px");
3443   proj[1] = hist->ProjectionY("vertex_check_py");
3444   proj[2] = hist->ProjectionZ("vertex_check_pz");
3445
3446   TFile::Open(data);
3447
3448   hist = (TH3*) gFile->Get("vertex_check");
3449
3450   TH1* proj2[3];
3451   proj2[0] = hist->ProjectionX("vertex_check_px2");
3452   proj2[1] = hist->ProjectionY("vertex_check_py2");
3453   proj2[2] = hist->ProjectionZ("vertex_check_pz2");
3454
3455   for (Int_t i=0; i<3; i++)
3456   {
3457     CompareQualityHists(proj[i], proj2[i], 1, 1);
3458   }
3459 }
3460
3461 void FitDiamondVsMult()
3462 {
3463   TFile::Open("analysis_esd_raw.root");
3464
3465   fVertexVsMult = (TH3*) gFile->Get("fVertexVsMult");
3466   fVertexVsMult->GetZaxis()->SetTitle("multiplicity");
3467
3468   TH2* proj[2];
3469   proj[0] = (TH2*) fVertexVsMult->Project3D("xz");
3470   proj[1] = (TH2*) fVertexVsMult->Project3D("yz");
3471
3472   gStyle->SetPadGridX(kTRUE);
3473   gStyle->SetPadGridY(kTRUE);
3474
3475   Int_t max = 40;
3476
3477   for (Int_t i=0; i<2; i++)
3478   {
3479     proj[i]->Rebin2D(4, 1);
3480     proj[i]->FitSlicesY();
3481
3482     c = new TCanvas(Form("c_%d", i), Form("c_%d", i), 800, 400);
3483     c->Divide(2, 1);
3484
3485     c->cd(1);
3486     hist = (TH1*) gROOT->FindObject(Form("fVertexVsMult_%sz_1", (i == 0) ? "x" : "y"));
3487     hist->GetXaxis()->SetRangeUser(0, max);
3488     hist->GetYaxis()->SetRangeUser(-0.4, 0.4);
3489     hist->Draw();
3490
3491     c->cd(2);
3492     hist = (TH1*) gROOT->FindObject(Form("fVertexVsMult_%sz_2", (i == 0) ? "x" : "y"));
3493     hist->GetXaxis()->SetRangeUser(0, max);
3494     hist->GetYaxis()->SetRangeUser(0, 0.2);
3495     hist->Draw();
3496
3497     c->SaveAs(Form("FitDiamondVsMult_%d.png", i));
3498   }
3499 }
3500
3501 void CompareQualityHists(const char* fileName1, const char* fileName2, const char* plotName, Int_t rebin1 = 1, Int_t rebin2 = 1, const char* exec = 0)
3502 {
3503   file1 = TFile::Open(fileName1);
3504   hist1 = (TH1*) file1->Get(plotName);
3505
3506   file2 = TFile::Open(fileName2);
3507   hist2 = (TH1*) file2->Get(plotName);
3508
3509   hist1->SetStats(0);
3510
3511   Printf("Entries in histograms: %f %f", hist1->Integral(), hist2->Integral());
3512
3513   if (exec)
3514   {
3515     hist1 = (TH1*) gROOT->ProcessLine(Form(exec, hist1, "hist1a"));
3516     hist2 = (TH1*) gROOT->ProcessLine(Form(exec, hist2, "hist2a"));
3517     hist1->Sumw2();
3518     hist2->Sumw2();
3519     Printf("Entries in histograms: %f %f", hist1->Integral(), hist2->Integral());
3520   }
3521
3522   CompareQualityHists(hist1, hist2, rebin1, rebin2);
3523 }
3524
3525 void CompareQualityHists(TH1* hist1, TH1* hist2, Int_t rebin1 = 1, Int_t rebin2 = 1)
3526 {
3527   hist1->SetLineColor(1);
3528   hist2->SetLineColor(2);
3529
3530   if (rebin1 != 0 && rebin2 != 0)
3531   {
3532     hist1->Rebin(TMath::Abs(rebin1));
3533     hist2->Rebin(TMath::Abs(rebin2));
3534   }
3535
3536   //hist2 = hist2->Rebin(hist1->GetNbinsX(), Form("%s_rebinned", hist2->GetName()), hist1->GetXaxis()->GetXbins()->GetArray());
3537
3538       //hist1->Scale(1.0 / 0.83);
3539
3540 //hist1->GetXaxis()->SetRangeUser(0, 50);
3541 /*  hist1->GetYaxis()->SetRangeUser(0.9, 1.2);
3542   hist1->Scale(1.0 / 0.808751);*/
3543
3544   //hist1->Scale(1.0 / 1.24632);
3545   //hist1->Scale(1.0 / 1.23821);
3546   //hist1->Scale(1.0 / 1.26213);
3547
3548   if (rebin1 > 0 && rebin2 > 0)
3549   {
3550     hist1->Scale(hist2->Integral() / hist1->Integral() * hist2->GetXaxis()->GetBinWidth(1) / hist1->GetXaxis()->GetBinWidth(1) / rebin1 * rebin2);
3551
3552     //hist1->Scale(0.5);
3553     //hist2->Scale(0.5);
3554   }
3555
3556   c = new TCanvas;
3557   if (strcmp(hist1->GetName(), "fDeltaTheta") == 0 || strcmp(hist1->GetName(), "fDeltaPhi") == 0 || strcmp(hist1->GetName(), "fMultVtx") == 0 || TString(hist1->GetName()).BeginsWith("vertex_check"))
3558     c->SetLogy();
3559
3560   if (TString(hist1->GetName()).BeginsWith("fMultiplicityESD"))
3561   {
3562     c->SetLogy();
3563     loadlibs();
3564     AliPWG0Helper::NormalizeToBinWidth(hist1);
3565     AliPWG0Helper::NormalizeToBinWidth(hist2);
3566   }
3567
3568   Printf("Means: %f %f %e", hist1->GetMean(), hist2->GetMean(), 1.0 - hist2->GetMean() / hist1->GetMean());
3569
3570   //hist1->GetYaxis()->SetRangeUser(0.01, hist1->GetMaximum() * 1.3);
3571   hist1->DrawCopy("HISTE");
3572   hist2->DrawCopy("HISTE SAME");
3573   gPad->SetGridx();
3574   gPad->SetGridy();
3575   //gPad->SetLogy();
3576   c->SaveAs(Form("%s_1.png", hist1->GetName()));
3577
3578   if (rebin1 == rebin2)
3579   {
3580     for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
3581       if (hist1->GetBinContent(i) == 0 && hist2->GetBinContent(i) > 0 || hist1->GetBinContent(i) > 0 && hist2->GetBinContent(i) == 0)
3582         Printf("Inconsistent bin %d: %f %f", i, hist1->GetBinContent(i), hist2->GetBinContent(i));
3583
3584     c2 = new TCanvas;
3585     hist1->GetYaxis()->SetRangeUser(0.5, 1.5);
3586     hist1->Divide(hist2);
3587     hist1->DrawCopy("HIST");
3588     gPad->SetGridx();
3589     gPad->SetGridy();
3590     c2->SaveAs(Form("%s_2.png", hist1->GetName()));
3591
3592     /*
3593     for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
3594       if (hist1->GetBinContent(i) > 0.9 && hist1->GetBinContent(i) < 1.1)
3595         hist1->SetBinContent(i, 0);
3596
3597     new TCanvas;
3598     hist1->SetMarkerStyle(20);
3599     hist1->DrawCopy("P");
3600     */
3601   }
3602 }
3603
3604 void DrawClustersVsTracklets()
3605 {
3606   TFile::Open("analysis_esd_raw.root");
3607
3608   hist = (TH2*) gFile->Get("fTrackletsVsClusters");
3609
3610   c = new TCanvas("c", "c", 600, 600);
3611   c->SetRightMargin(0.05);
3612   c->SetTopMargin(0.05);
3613
3614   hist->SetStats(0);
3615   hist->GetYaxis()->SetRangeUser(0, 400);
3616   hist->GetYaxis()->SetTitleOffset(1.3);
3617   hist->GetXaxis()->SetRangeUser(0, 30);
3618   hist->Draw("BOX");
3619
3620   func = new TF1("func", "80 + x * 11", 0, 30);
3621   func->Draw("SAME");
3622
3623   c->SaveAs("clusters_vs_tracklets.eps");
3624 }
3625
3626 void VertexPlotBackgroundNote()
3627 {
3628   TFile::Open("all.root");
3629
3630   hist = (TH3*) gFile->Get("vertex_check");
3631   proj = (TH1*) hist->ProjectionZ()->Clone("all");
3632   proj->Rebin(2);
3633
3634   proj->Draw();
3635
3636   TFile::Open("analysis_esd_raw.root");
3637   hist = (TH3*) gFile->Get("vertex_check");
3638   proj = (TH1*) hist->ProjectionZ()->Clone("afterbg");
3639   proj->Rebin(2);
3640
3641   proj->SetLineColor(2);
3642   proj->Draw("SAME");
3643 }
3644
3645 void BackgroundAnalysis(const char* signal, const char* background)
3646 {
3647   TFile::Open(signal);
3648   signalHist = (TH2*) gFile->Get("fTrackletsVsClusters");
3649
3650   TFile::Open(background);
3651   backgroundHist = (TH2*) gFile->Get("fTrackletsVsClusters");
3652
3653   Printf("For events with >= 1 tracklet:");
3654
3655   func = new TF1("func", "[0] + x * 11", 0, 30);
3656   for (Int_t a = 50; a <= 100; a += 10)
3657   {
3658     func->SetParameter(0, a);
3659
3660     Float_t signalCount = 0;
3661     Float_t backgroundCount = 0;
3662     for (Int_t x = 2; x <= signalHist->GetNbinsX(); x++)
3663     {
3664       signalCount += signalHist->Integral(x, x, signalHist->GetYaxis()->FindBin(func->Eval(signalHist->GetXaxis()->GetBinCenter(x))), signalHist->GetNbinsY());
3665       backgroundCount += backgroundHist->Integral(x, x, signalHist->GetYaxis()->FindBin(func->Eval(signalHist->GetXaxis()->GetBinCenter(x))), signalHist->GetNbinsY());
3666     }
3667
3668     Float_t signalFraction = 100.0 * signalCount / signalHist->Integral(2, signalHist->GetNbinsX(), 1, signalHist->GetNbinsY());
3669     Float_t backgroundFraction = 100.0 * backgroundCount / backgroundHist->Integral(2, signalHist->GetNbinsX(), 1, signalHist->GetNbinsY());
3670
3671     Printf("Cut at a = %d; Removed %.2f %% of the background (%.0f events); Removed %.2f %% of the signal", a, backgroundFraction, backgroundCount, signalFraction);
3672   }
3673 }
3674
3675 void ZPhiPlots()
3676 {
3677   TFile::Open("analysis_esd_raw.root");
3678
3679   for (Int_t i=0; i<2; i++)
3680   {
3681     hist = (TH2*) gFile->Get(Form("fZPhi_%d", i));
3682
3683     c = new TCanvas;
3684     hist->SetStats(0);
3685     hist->Draw("COLZ");
3686     c->SaveAs(Form("ZPhi_%d.png", i));
3687   }
3688 }
3689
3690 void DrawStats(Bool_t all = kFALSE)
3691 {
3692   if (all)
3693   {
3694     Int_t count = 4;
3695     const char* list[] = { "CINT1B-ABCE-NOPF-ALL/spd", "CINT1A-ABCE-NOPF-ALL/spd", "CINT1C-ABCE-NOPF-ALL/spd", "CINT1-E-NOPF-ALL/spd" };
3696   }
3697   else
3698   {
3699     Int_t count = 1;
3700     const char* list[] = { "." };
3701   }
3702
3703   for (Int_t i=0; i<count; i++)
3704   {
3705     TFile::Open(Form("%s/analysis_esd_raw.root", list[i]));
3706
3707     hist = (TH2*) gFile->Get("fStats2");
3708
3709     c = new TCanvas(list[i], list[i], 800, 600);
3710     gPad->SetBottomMargin(0.2);
3711     gPad->SetLeftMargin(0.2);
3712     gPad->SetRightMargin(0.2);
3713     hist->Draw("TEXT");
3714     hist->SetMarkerSize(2);
3715     //hist->GetYaxis()->SetRangeUser(0, 0);
3716
3717     gROOT->Macro("increaseFonts.C");
3718
3719     c->SaveAs(Form("%s/stats.png", list[i]));
3720   }
3721 }
3722
3723 void CompareMCDataTrigger(const char* mcFile, const char* dataFile)
3724 {
3725   TH1* stat[2];
3726
3727   TFile::Open(mcFile);
3728   mc = (TH1*) gFile->Get("trigger_histograms_/fHistFiredBitsSPD");
3729   stat[0] = (TH1*) gFile->Get("fHistStatistics");
3730
3731   TFile::Open(dataFile);
3732   data = (TH1*) gFile->Get("trigger_histograms_+CINT1B-ABCE-NOPF-ALL/fHistFiredBitsSPD");
3733   if (!data)
3734     data = (TH1*) gFile->Get("trigger_histograms_+CSMBB-ABCE-NOPF-ALL/fHistFiredBitsSPD");
3735
3736   stat[1] = (TH1*) gFile->Get("fHistStatistics");
3737
3738   CompareQualityHists(mc, data);
3739
3740   for (Int_t i=0; i<2; i++)
3741   {
3742     Float_t total = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("Trigger class"), 1);
3743     Float_t spd = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("FO >= 2"), 1);
3744     Float_t v0A = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("V0A"), 1);
3745     Float_t v0C = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("V0C"), 1);
3746
3747     Printf("%s:\nSPD / V0A: %.3f\nSPD / V0C: %.3f\nV0A / V0C: %.3f", (i == 0) ? "MC  " : "Data", spd / v0A, spd / v0C, v0A / v0C);
3748     Printf("SPD / Total: %.3f\nV0A / Total: %.3f\nV0C / Total: %.3f\n", spd / total, v0A / total, v0C / total);
3749   }
3750 }
3751
3752 void CompareMCDatadNdEta(const char* mcFile, const char* dataFile)
3753 {
3754   //CompareQualityHists(mcFile, dataFile, "fEtaPhi", 4, 4, "((TH2*)%p)->ProjectionY(\"%s\", 1, 40)");
3755   //CompareQualityHists(mcFile, dataFile, "fEtaPhi", 4, 4, "((TH2*)%p)->ProjectionY(\"%s\", 41, 80)");
3756
3757   CompareQualityHists(mcFile, dataFile, "fEtaPhi", 1, 1, "((TH2*)%p)->ProjectionX(\"%s\", 271, 360)");
3758 }
3759
3760 void TrigVsTrigVtx(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
3761 {
3762   loadlibs();
3763   if (!TFile::Open(fileName))
3764     return;
3765
3766   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
3767   if (!dNdEtaCorrection->LoadHistograms())
3768     return;
3769
3770   TH2* eTrig =    dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
3771   TH2* eTrigVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
3772
3773   eTrig_proj = eTrig->ProjectionY("y1", eTrig->GetYaxis()->FindBin(-9.9), eTrig->GetYaxis()->FindBin(9.9));
3774   eTrigVtx_proj = eTrigVtx->ProjectionY("y2", eTrig->GetYaxis()->FindBin(-9.9), eTrig->GetYaxis()->FindBin(9.9));
3775
3776   new TCanvas;
3777   eTrig_proj->Draw();
3778   eTrig_proj->GetXaxis()->SetRangeUser(0, 20);
3779   eTrigVtx_proj->SetLineColor(2);
3780   eTrigVtx_proj->Draw("SAME");
3781
3782   gPad->SetLogy();
3783 }
3784
3785 void PrintAverageNSDCorrectionFactors()
3786 {
3787   // factors estimated from MC, can be slighly different with data b/c correction is applies as function of measured multiplicity
3788
3789   loadlibs();
3790
3791   if (!TFile::Open("correction_mapprocess-types.root"))
3792     return;
3793
3794   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND");
3795   if (!dNdEtaCorrection->LoadHistograms())
3796     return;
3797
3798   AlidNdEtaCorrection* dNdEtaCorrection2 = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD");
3799   if (!dNdEtaCorrection2->LoadHistograms())
3800     return;
3801
3802   // for scaling factors see drawSystematics.C; GetRelativeFractions()
3803   // 900 GeV
3804   //dNdEtaCorrection->Scale(1.06);
3805   //dNdEtaCorrection->Add(dNdEtaCorrection2, 9.5 / 12.3);
3806   // 2.36 TeV
3807   dNdEtaCorrection->Scale(1.036);
3808   dNdEtaCorrection->Add(dNdEtaCorrection2, 0.075 * 1.43 / 0.127);
3809
3810   Printf("event adding: %f", dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral());
3811
3812   Printf("track adding: %f", dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetMeasuredHistogram()->Integral());
3813
3814   AlidNdEtaCorrection* dNdEtaCorrection3 = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD");
3815   if (!dNdEtaCorrection3->LoadHistograms())
3816     return;
3817
3818   // 900 GeV
3819   //dNdEtaCorrection3->Scale(0.153 / 0.189);
3820   // 2.36 TeV
3821   dNdEtaCorrection3->Scale(0.159 / 0.166);
3822   dNdEtaCorrection->Add(dNdEtaCorrection3);
3823
3824   Printf("event subtraction: %f", dNdEtaCorrection3->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral());
3825
3826   Printf("track subtraction: %f", dNdEtaCorrection3->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetMeasuredHistogram()->Integral());
3827
3828   dNdEtaCorrection->GetTriggerBiasCorrectionNSD()->PrintInfo(0.0);
3829 }
3830