PWG0/dNdEta/drawPlots.C

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