1590528538f587e155855bdd7d0a184ffd6a1680
[u/mrichter/AliRoot.git] / 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 TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
678 {
679   gSystem->Load("libPWG0base");
680
681   TFile::Open(fileName);
682   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
683   dNdEtaCorrection->LoadHistograms();
684
685   TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
686   TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
687
688   Prepare2DPlot(corrTrigger);
689   corrTrigger->SetTitle("b) Trigger bias correction");
690
691   Prepare2DPlot(corrVtx);
692   corrVtx->SetTitle("a) Vertex reconstruction correction");
693
694   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
695   corrVtx->GetYaxis()->SetTitle("Multiplicity");
696
697   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
698   canvas->Divide(2, 1);
699
700   canvas->cd(1);
701   InitPadCOLZ();
702   corrVtx->DrawCopy("COLZ");
703
704   canvas->cd(2);
705   InitPadCOLZ();
706   corrTrigger->DrawCopy("COLZ");
707
708   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
709   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
710
711   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
712   canvas->Divide(2, 1);
713
714   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
715   corrVtx->GetYaxis()->SetRangeUser(0, 5);
716
717   canvas->cd(1);
718   InitPadCOLZ();
719   corrVtx->DrawCopy("COLZ");
720
721   canvas->cd(2);
722   InitPadCOLZ();
723   corrTrigger->DrawCopy("COLZ");
724
725   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
726   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
727 }
728
729 void TriggerBias(const char* fileName = "correction_map.root")
730 {
731   TFile* file = TFile::Open(fileName);
732
733   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
734
735   Prepare2DPlot(corr);
736   corr->SetTitle("Trigger bias correction");
737
738   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
739   InitPadCOLZ();
740   corr->DrawCopy("COLZ");
741
742   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
743   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
744
745   corr->GetYaxis()->SetRangeUser(0, 5);
746
747   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
748   InitPadCOLZ();
749   corr->DrawCopy("COLZ");
750
751   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
752   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
753 }
754
755 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
756 {
757   gSystem->Load("libPWG0base");
758
759   TFile* file = TFile::Open(fileName);
760   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
761   dNdEtaCorrection->LoadHistograms();
762
763   TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
764   TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
765
766   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
767   canvas->Divide(2, 1);
768
769   canvas->cd(1);
770   InitPad();
771
772   Prepare1DPlot(hist);
773   hist->SetTitle("");
774   hist->GetYaxis()->SetTitle("correction factor");
775   hist->GetYaxis()->SetRangeUser(1, 1.5);
776   hist->GetYaxis()->SetTitleOffset(1.6);
777   hist->Draw();
778
779   canvas->cd(2);
780   InitPad();
781
782   Prepare1DPlot(hist2);
783   hist2->SetTitle("");
784   hist2->GetYaxis()->SetTitle("correction factor");
785   hist2->GetXaxis()->SetRangeUser(0, 5);
786   hist2->GetYaxis()->SetTitleOffset(1.6);
787   hist2->GetXaxis()->SetTitle("multiplicity");
788   hist2->Draw();
789
790   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
791   pave->SetFillColor(0);
792   pave->AddText("|z| < 10 cm");
793   pave->Draw();
794
795   canvas->SaveAs("TriggerBias1D.eps");
796 }
797
798 void VtxRecon()
799 {
800   TFile* file = TFile::Open("correction_map.root");
801
802   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
803
804   Prepare2DPlot(corr);
805   corr->SetTitle("Vertex reconstruction correction");
806
807   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
808   InitPadCOLZ();
809   corr->DrawCopy("COLZ");
810
811   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
812   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
813
814   corr->GetYaxis()->SetRangeUser(0, 5);
815
816   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
817   InitPadCOLZ();
818   corr->DrawCopy("COLZ");
819
820   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
821   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
822 }
823
824 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
825 {
826   gSystem->Load("libPWG0base");
827
828   TFile* file = TFile::Open(fileName);
829   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
830   dNdEtaCorrection->LoadHistograms();
831
832   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
833   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
834
835   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
836   canvas->Divide(2, 1);
837
838   canvas->cd(1);
839   InitPad();
840
841   Prepare1DPlot(hist);
842   hist->SetTitle("");
843   hist->GetYaxis()->SetTitle("correction factor");
844   hist->GetYaxis()->SetRangeUser(1, 1.8);
845   hist->GetYaxis()->SetTitleOffset(1.6);
846   hist->DrawCopy();
847
848   canvas->cd(2);
849   InitPad();
850
851   Prepare1DPlot(hist2);
852   hist2->SetTitle("");
853   hist2->GetYaxis()->SetTitle("correction factor");
854   hist2->GetXaxis()->SetRangeUser(0, 20);
855   hist2->GetYaxis()->SetTitleOffset(1.6);
856   hist2->GetXaxis()->SetTitle("multiplicity");
857   hist2->Draw();
858
859   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
860   pave->SetFillColor(0);
861   pave->AddText("|z| < 10 cm");
862   pave->Draw();
863
864   canvas->SaveAs("VtxRecon1D.eps");
865
866   Correction1DCreatePlots(fileName, folderName, 9.9, 2);
867
868   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
869   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
870
871   Prepare1DPlot(corrX);
872   Prepare1DPlot(corrZ);
873
874   corrX->GetYaxis()->SetTitleOffset(1.5);
875   corrZ->GetYaxis()->SetTitleOffset(1.5);
876
877   corrX->SetTitle("a) z projection");
878   corrZ->SetTitle("b) p_{T} projection");
879
880   corrX->GetYaxis()->SetTitle("correction factor");
881   corrZ->GetYaxis()->SetTitle("correction factor");
882
883   corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
884
885   TString canvasName;
886   canvasName.Form("VtxRecon1D_Track");
887   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
888   canvas->Divide(2, 1);
889
890   canvas->cd(1);
891   InitPad();
892   corrX->DrawCopy();
893
894   canvas->cd(2);
895   InitPad();
896   gPad->SetLogx();
897   corrZ->Draw();
898
899   canvas->SaveAs("VtxRecon1D_Track.eps");
900   canvas->SaveAs("VtxRecon1D_Track.gif");
901 }
902
903 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
904 {
905   gSystem->Load("libPWG0base");
906
907   TFile::Open(fileName);
908   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
909   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
910
911   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
912   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
913
914   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
915   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
916
917   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
918   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
919   gene->GetXaxis()->SetRangeUser(-10, 10);
920   meas->GetXaxis()->SetRangeUser(-10, 10);
921
922   Float_t eff1 = gene->Integral() / meas->Integral();
923   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
924
925   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
926
927   gene->GetZaxis()->SetRangeUser(0.3, 10);
928   meas->GetZaxis()->SetRangeUser(0.3, 10);
929
930   Float_t eff2 = gene->Integral() / meas->Integral();
931   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
932
933   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
934
935   gene->GetZaxis()->SetRangeUser(0.3, 1);
936   meas->GetZaxis()->SetRangeUser(0.3, 1);
937
938   Float_t eff3 = gene->Integral() / meas->Integral();
939   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
940
941   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
942 }
943
944 void Correction1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType = 0)
945 {
946   TFile::Open(fileName);
947   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
948   dNdEtaCorrection->LoadHistograms();
949
950   TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
951   TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetMeasuredHistogram();
952
953   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
954   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
955   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
956   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
957   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kFALSE);
958   gene->GetYaxis()->SetRange(0, 0);
959   meas->GetYaxis()->SetRange(0, 0);
960
961   gene->GetXaxis()->SetRangeUser(-9.9, 9.9);
962   meas->GetXaxis()->SetRangeUser(-9.9, 9.9);
963   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kFALSE);
964   gene->GetZaxis()->SetRange(0, 0);
965   meas->GetZaxis()->SetRange(0, 0);
966
967   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
968   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
969   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kFALSE);
970 }
971
972 TCanvas* Correction1D(Int_t correctionType = 0, const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
973 {
974   gSystem->Load("libPWG0base");
975
976   Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType);
977
978   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
979   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
980   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
981
982   Prepare1DPlot(corrX);
983   Prepare1DPlot(corrY);
984   Prepare1DPlot(corrZ);
985
986   corrX->SetTitle("a) z projection");
987   corrY->SetTitle("b) #eta projection");
988   corrZ->SetTitle("c) p_{T} projection");
989
990   corrX->GetYaxis()->SetTitle("correction factor");
991   corrY->GetYaxis()->SetTitle("correction factor");
992   corrZ->GetYaxis()->SetTitle("correction factor");
993   corrX->GetYaxis()->SetTitleOffset(1.7);
994   corrY->GetYaxis()->SetTitleOffset(1.7);
995   corrZ->GetYaxis()->SetTitleOffset(1.7);
996   corrX->GetYaxis()->SetRangeUser(0.8, 1.5);
997   corrY->GetYaxis()->SetRangeUser(0.8, 1.5);
998   corrZ->GetYaxis()->SetRangeUser(0.8, 1.5);
999
1000   corrZ->GetXaxis()->SetRangeUser(0.11, upperPtLimit);
1001
1002   TString canvasName;
1003   canvasName.Form(Form("Correction1D_%d_%s_%f", correctionType, fileName, upperPtLimit));
1004   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1005   canvas->Divide(3, 1);
1006
1007   TLatex* Tl = new TLatex;
1008   Tl->SetTextSize(0.04);
1009   Tl->SetBit(TLatex::kTextNDC);
1010
1011   canvas->cd(1);
1012   InitPad();
1013   corrX->DrawCopy();
1014   Tl->DrawLatex(0.6, 0.8, "0.3 < p_{T} < 10");
1015   Tl->DrawLatex(0.6, 0.75, "|#eta| < 0.8");
1016
1017   canvas->cd(2);
1018   InitPad();
1019   corrY->Draw();
1020   Tl->DrawLatex(0.6, 0.8, "0.3 < p_{T} < 10");
1021   Tl->DrawLatex(0.6, 0.75, "|vtx-z| < 10");
1022
1023   canvas->cd(3);
1024   InitPad();
1025   gPad->SetLogx();
1026   corrZ->Draw();
1027   Tl->DrawLatex(0.6, 0.8, "|vtx-z| < 10");
1028   Tl->DrawLatex(0.6, 0.75, "|#eta| < 0.8");
1029
1030   return canvas;
1031 }
1032
1033 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1034 {
1035   gSystem->Load("libPWG0base");
1036
1037   Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
1038
1039   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1040   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1041   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1042
1043   Prepare1DPlot(corrX);
1044   Prepare1DPlot(corrY);
1045   Prepare1DPlot(corrZ);
1046
1047   corrX->SetTitle("a) z projection");
1048   corrY->SetTitle("a) #eta projection");
1049   corrZ->SetTitle("b) p_{T} projection");
1050
1051   corrY->GetYaxis()->SetTitle("correction factor");
1052   corrZ->GetYaxis()->SetTitle("correction factor");
1053
1054   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1055
1056   TString canvasName;
1057   canvasName.Form("Track2Particle1D_%s", folder);
1058   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1059   canvas->Divide(3, 1);
1060
1061   canvas->cd(1);
1062   InitPad();
1063   corrX->DrawCopy();
1064
1065   canvas->cd(2);
1066   InitPad();
1067   corrY->Draw();
1068
1069   canvas->cd(3);
1070   InitPad();
1071   corrZ->Draw();
1072
1073   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
1074   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
1075
1076   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
1077
1078   canvasName.Form("Track2Particle1D_%s_etapt", folder);
1079   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
1080   canvas->Divide(2, 1);
1081
1082   canvas->cd(1);
1083   InitPad();
1084   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
1085   corrY->GetYaxis()->SetRangeUser(1, 1.5);
1086   corrY->GetYaxis()->SetTitleOffset(1.5);
1087   corrY->DrawCopy();
1088   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
1089   pave->AddText("|z| < 10 cm");
1090   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
1091   pave->Draw();
1092
1093   canvas->cd(2);
1094   InitPad();
1095   gPad->SetLogx();
1096   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
1097   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
1098   corrZ->GetYaxis()->SetTitleOffset(1.5);
1099   corrZ->DrawCopy();
1100   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
1101   pave->AddText("|z| < 10 cm");
1102   pave->AddText("|#eta| < 0.8");
1103   pave->Draw();
1104
1105   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
1106   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
1107 }
1108
1109 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1110 {
1111   gSystem->Load("libPWG0base");
1112
1113   // particle type
1114   for (Int_t particle=0; particle<4; ++particle)
1115   {
1116     TString dirName;
1117     dirName.Form("correction_%d", particle);
1118     Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
1119
1120     TString tmpx, tmpy, tmpz;
1121     tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1122     tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1123     tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
1124
1125     TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1126     TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1127     TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
1128
1129     Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
1130
1131     TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1132     TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1133     TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
1134
1135     posX->Divide(negX);
1136     posY->Divide(negY);
1137     posZ->Divide(negZ);
1138
1139     Prepare1DPlot(posX);
1140     Prepare1DPlot(posY);
1141     Prepare1DPlot(posZ);
1142
1143     Float_t min = 0.8;
1144     Float_t max = 1.2;
1145
1146     posX->SetMinimum(min);
1147     posX->SetMaximum(max);
1148     posY->SetMinimum(min);
1149     posY->SetMaximum(max);
1150     posZ->SetMinimum(min);
1151     posZ->SetMaximum(max);
1152
1153     posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1154
1155     posX->GetYaxis()->SetTitleOffset(1.7);
1156     posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1157     posY->GetYaxis()->SetTitleOffset(1.7);
1158     posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1159     posZ->GetYaxis()->SetTitleOffset(1.7);
1160     posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
1161
1162     posZ->GetXaxis()->SetRangeUser(0, 1);
1163
1164     TString canvasName;
1165     canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
1166
1167     TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1168     canvas->Divide(3, 1);
1169
1170     canvas->cd(1);
1171     InitPad();
1172     posX->DrawCopy();
1173
1174     canvas->cd(2);
1175     InitPad();
1176     posY->DrawCopy();
1177
1178     canvas->cd(3);
1179     InitPad();
1180     posZ->DrawCopy();
1181
1182     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1183     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1184   }
1185 }
1186
1187 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
1188 {
1189   TFile::Open(fileName);
1190   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1191   dNdEtaCorrection->LoadHistograms();
1192
1193   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1194   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1195
1196   gene->GetZaxis()->SetRangeUser(0.3, 10);
1197   meas->GetZaxis()->SetRangeUser(0.3, 10);
1198   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
1199   gene->GetZaxis()->SetRange(0, 0);
1200   meas->GetZaxis()->SetRange(0, 0);
1201
1202   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1203   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1204   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
1205   gene->GetYaxis()->SetRange(0, 0);
1206   meas->GetYaxis()->SetRange(0, 0);
1207
1208   gene->GetXaxis()->SetRangeUser(-10, 10);
1209   meas->GetXaxis()->SetRangeUser(-10, 10);
1210   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
1211   gene->GetXaxis()->SetRange(0, 0);
1212   meas->GetXaxis()->SetRange(0, 0);
1213 }
1214
1215 TCanvas* Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1216 {
1217   gSystem->Load("libPWG0base");
1218
1219   Track2Particle2DCreatePlots(fileName);
1220
1221   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1222   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1223   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
1224
1225   Prepare2DPlot(corrYX);
1226   Prepare2DPlot(corrZX);
1227   Prepare2DPlot(corrZY);
1228
1229   const char* title = "";
1230   corrYX->SetTitle(title);
1231   corrZX->SetTitle(title);
1232   corrZY->SetTitle(title);
1233
1234   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1235   canvas->Divide(3, 1);
1236
1237   canvas->cd(1);
1238   InitPadCOLZ();
1239   corrYX->Draw("COLZ");
1240
1241   canvas->cd(2);
1242   InitPadCOLZ();
1243   corrZX->Draw("COLZ");
1244
1245   canvas->cd(3);
1246   InitPadCOLZ();
1247   corrZY->Draw("COLZ");
1248
1249   canvas->SaveAs(Form("spd_corr_track2particle_%d.gif", gMax));
1250   canvas->SaveAs(Form("spd_corr_track2particle_%d.eps", gMax));
1251
1252   return canvas;
1253 }
1254
1255 void CompareTrack2Particle2D()
1256 {
1257   gSystem->Load("libPWG0base");
1258
1259   Track2Particle2DCreatePlots("correction_maponly-positive.root");
1260
1261   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("pos_yx"));
1262   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("pos_zx"));
1263   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("pos_zy"));
1264
1265   Track2Particle2DCreatePlots("correction_maponly-negative.root");
1266
1267   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("neg_yx"));
1268   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("neg_zx"));
1269   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("neg_zy"));
1270
1271   posYX->Divide(negYX);
1272   posZX->Divide(negZX);
1273   posZY->Divide(negZY);
1274
1275   Prepare2DPlot(posYX);
1276   Prepare2DPlot(posZX);
1277   Prepare2DPlot(posZY);
1278
1279   Float_t min = 0.8;
1280   Float_t max = 1.2;
1281
1282   posYX->SetMinimum(min);
1283   posYX->SetMaximum(max);
1284   posZX->SetMinimum(min);
1285   posZX->SetMaximum(max);
1286   posZY->SetMinimum(min);
1287   posZY->SetMaximum(max);
1288
1289   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1290   canvas->Divide(3, 1);
1291
1292   canvas->cd(1);
1293   InitPadCOLZ();
1294   posYX->Draw("COLZ");
1295
1296   canvas->cd(2);
1297   InitPadCOLZ();
1298   posZX->Draw("COLZ");
1299
1300   canvas->cd(3);
1301   InitPadCOLZ();
1302   posZY->Draw("COLZ");
1303
1304   canvas->SaveAs("CompareTrack2Particle2D.gif");
1305   canvas->SaveAs("CompareTrack2Particle2D.eps");
1306 }
1307
1308 void Track2Particle3D()
1309 {
1310   // get left margin proper
1311
1312   TFile* file = TFile::Open("correction_map.root");
1313
1314   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1315
1316   corr->SetTitle("Correction Factor");
1317   SetRanges(corr->GetZaxis());
1318
1319   Prepare3DPlot(corr);
1320
1321   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1322   canvas->SetTheta(29.428);
1323   canvas->SetPhi(16.5726);
1324
1325   corr->Draw();
1326
1327   canvas->SaveAs("Track2Particle3D.gif");
1328   canvas->SaveAs("Track2Particle3D.eps");
1329 }
1330
1331 void Track2Particle3DAll()
1332 {
1333   TFile* file = TFile::Open("correction_map.root");
1334
1335   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1336   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1337   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1338
1339   gene->SetTitle("Generated Particles");
1340   meas->SetTitle("Measured Tracks");
1341   corr->SetTitle("Correction Factor");
1342
1343   Prepare3DPlot(gene);
1344   Prepare3DPlot(meas);
1345   Prepare3DPlot(corr);
1346
1347   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1348   canvas->Divide(3, 1);
1349
1350   canvas->cd(1);
1351   InitPad();
1352   gene->Draw();
1353
1354   canvas->cd(2);
1355   meas->Draw();
1356
1357   canvas->cd(3);
1358   corr->Draw();
1359
1360   canvas->SaveAs("Track2Particle3DAll.gif");
1361   canvas->SaveAs("Track2Particle3DAll.eps");
1362 }
1363
1364 void MultiplicityMC(Int_t xRangeMax = 50)
1365 {
1366   TFile* file = TFile::Open("multiplicityMC.root");
1367
1368   if (!file)
1369   {
1370     printf("multiplicityMC.root could not be opened.\n");
1371     return;
1372   }
1373
1374   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1375   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1376   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1377
1378   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1379   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1380   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1381   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1382   {
1383     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1384     proj->Fit("gaus", "0");
1385     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1386     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1387
1388     continue;
1389
1390     // draw for debugging
1391     new TCanvas;
1392     proj->DrawCopy();
1393     proj->GetFunction("gaus")->DrawCopy("SAME");
1394   }
1395
1396   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1397
1398   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1399   {
1400     Float_t mean = correction->GetBinContent(i);
1401     Float_t width = correctionWidth->GetBinContent(i);
1402
1403     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1404     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1405     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1406
1407     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1408     {
1409       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1410     }
1411   }
1412
1413   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1414   fMultiplicityESDCorrectedRebinned->Rebin(10);
1415   fMultiplicityESDCorrectedRebinned->Scale(0.1);
1416
1417   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1418   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1419   ratio->Divide(fMultiplicityMC);
1420
1421   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1422   ratio2->Divide(fMultiplicityMC);
1423
1424   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1425   canvas->Divide(3, 2);
1426
1427   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1428   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1429   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1430   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1431   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1432   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1433   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1434
1435   canvas->cd(1); //InitPad();
1436   fMultiplicityESD->Draw();
1437   fMultiplicityMC->SetLineColor(2);
1438   fMultiplicityMC->Draw("SAME");
1439
1440   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1441   legend->AddEntry(fMultiplicityESD, "ESD");
1442   legend->AddEntry(fMultiplicityMC, "MC");
1443   legend->Draw();
1444
1445   canvas->cd(2);
1446   fCorrelation->Draw("COLZ");
1447
1448   canvas->cd(3);
1449   correction->Draw();
1450   //correction->Fit("pol1");
1451   correctionWidth->SetLineColor(2);
1452   correctionWidth->Draw("SAME");
1453
1454   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1455   legend->AddEntry(correction, "#bar{x}");
1456   legend->AddEntry(correctionWidth, "#sigma");
1457   legend->Draw();
1458
1459   canvas->cd(4);
1460   ratio->Draw();
1461
1462   ratio2->SetLineColor(2);
1463   ratio2->Draw("SAME");
1464
1465   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1466   legend->AddEntry(ratio, "uncorrected");
1467   legend->AddEntry(ratio2, "corrected");
1468   legend->Draw();
1469
1470   canvas->cd(5);
1471   fMultiplicityESDCorrected->SetLineColor(kBlue);
1472   fMultiplicityESDCorrected->Draw();
1473   fMultiplicityMC->Draw("SAME");
1474   fMultiplicityESD->Draw("SAME");
1475
1476   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1477   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1478   legend->AddEntry(fMultiplicityMC, "MC");
1479   legend->AddEntry(fMultiplicityESD, "ESD");
1480   legend->Draw();
1481
1482   canvas->cd(6);
1483   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
1484   fMultiplicityESDCorrectedRebinned->Draw();
1485   fMultiplicityMC->Draw("SAME");
1486
1487   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1488   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1489   legend->AddEntry(fMultiplicityMC, "MC");
1490   legend->Draw();
1491
1492   canvas->SaveAs("MultiplicityMC.gif");
1493 }
1494
1495 void MultiplicityESD()
1496 {
1497   TFile* file = TFile::Open("multiplicityESD.root");
1498
1499   if (!file)
1500   {
1501     printf("multiplicityESD.root could not be opened.\n");
1502     return;
1503   }
1504
1505   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1506
1507   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1508
1509   fMultiplicityESD->Draw();
1510 }
1511
1512 void drawPlots(Int_t max)
1513 {
1514   gMax = max;
1515
1516   ptCutoff();
1517   TriggerBias();
1518   VtxRecon();
1519   Track2Particle2D();
1520   Track2Particle3D();
1521   ptSpectrum();
1522   dNdEta();
1523 }
1524
1525 void drawPlots()
1526 {
1527   drawPlots(5);
1528   drawPlots(2);
1529 }
1530
1531 void CompareCorrection2Measured(const char* dataInput = "analysis_esd_raw.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1532 {
1533   loadlibs();
1534
1535   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1536   TFile::Open(correctionMapFile);
1537   dNdEtaCorrection->LoadHistograms();
1538
1539   TFile* file = TFile::Open(dataInput);
1540
1541   if (!file)
1542   {
1543     cout << "Error. File not found" << endl;
1544     return;
1545   }
1546
1547   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1548   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1549
1550   gROOT->cd();
1551   
1552   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
1553   hist1->SetTitle("mc");
1554   Printf("mc contains %f entries", hist1->Integral());
1555   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()));
1556
1557   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
1558   hist2->SetTitle("esd");
1559   Printf("esd contains %f entries", hist2->Integral());
1560   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()));
1561
1562   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1563   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
1564
1565   hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1566   hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1567   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
1568
1569   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
1570   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
1571   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
1572   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
1573   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
1574 }
1575
1576 void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
1577 {
1578   loadlibs();
1579
1580   TFile* file = TFile::Open(dataInput);
1581
1582   if (!file)
1583   {
1584     cout << "Error. File not found" << endl;
1585     return;
1586   }
1587
1588   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1589   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1590
1591   TFile* file = TFile::Open(dataInput2);
1592
1593   if (!file)
1594   {
1595     cout << "Error. File not found" << endl;
1596     return;
1597   }
1598
1599   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
1600   fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
1601
1602   gROOT->cd();
1603
1604   TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
1605   hist1->SetTitle("esd1");
1606   Printf("esd1 contains %f entries", hist1->GetEntries());
1607   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()));
1608
1609   TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
1610   hist2->SetTitle("esd2");
1611   Printf("esd2 contains %f entries", hist2->GetEntries());
1612   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()));
1613
1614   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1615
1616   new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
1617   new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
1618   new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
1619
1620   TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
1621   TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
1622
1623   Printf("event1 contains %f entries", event1->GetEntries());
1624   Printf("event2 contains %f entries", event2->GetEntries());
1625   Printf("event1 integral is %f", event1->Integral());
1626   Printf("event2 integral is %f", event2->Integral());
1627   Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
1628   Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
1629
1630   projx1 = event1->ProjectionX();
1631   projx2 = event2->ProjectionX();
1632
1633   new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
1634
1635   projx1->Divide(projx2);
1636   new TCanvas; projx1->Draw();
1637
1638   event1->Divide(event2);
1639   new TCanvas; event1->Draw("COLZ");
1640
1641 }
1642
1643 void DrawTrackletOrigin()
1644 {
1645   TFile::Open("correction_map.root");
1646
1647   Int_t colors[]  = {1,2,3,4,6,7,8,102};
1648
1649   Int_t maxHists = 8;
1650   TH1* hist[8];
1651
1652   const char* titles[] = { "PP", "SS", "PP'", "PS", "PS*", "SP", "SS'", "" };
1653
1654   TLegend* legend = new TLegend(0.75, 0.6, 0.95, 0.95);
1655
1656   Int_t total = 0;
1657   for (Int_t i=0; i<maxHists; i++)
1658   {
1659     hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
1660     //hist[i]->Rebin(20);
1661     hist[i]->SetStats(kFALSE);
1662     hist[i]->SetLineColor(colors[i]);
1663     hist[i]->GetXaxis()->SetRangeUser(-0.2, 0.2);
1664     hist[i]->Draw(((i == 0) ? "" : "SAME"));
1665
1666     total += hist[i]->GetEntries();
1667
1668     if (i != 7)
1669       legend->AddEntry(hist[i], titles[i]);
1670   }
1671
1672   legend->Draw();
1673   gPad->SetLogy();
1674
1675   Printf("Total: %d", total);
1676   for (Int_t i=0; i<maxHists; i++)
1677     Printf("Histogram %d (%s) containts %.2f %% of the entries", i, titles[i], 100.0 * hist[i]->GetEntries() / total);
1678
1679   printf("|  Delta phi  |  Acc. %%  |  ");
1680   for (Int_t i=0; i<maxHists; i++)
1681     printf("%3s %%   |  ", titles[i]);
1682   Printf("");
1683
1684   for (Float_t f = 0.01; f < 0.09; f += 0.01)
1685   {
1686     Int_t integralBegin = hist[0]->GetXaxis()->FindBin(-f);
1687     Int_t integralEnd = hist[0]->GetXaxis()->FindBin(f);
1688
1689     Int_t total2 = 0;
1690     for (Int_t i=0; i<maxHists; i++)
1691       total2 += (Int_t) hist[i]->Integral(integralBegin, integralEnd);
1692
1693     printf("|    %.2f     |  %6.2f  |  ", f, 100.0 * total2 / total);
1694
1695     for (Int_t i=0; i<maxHists; i++)
1696       printf("%6.2f  |  ", (hist[i]->GetEntries() > 0) ? (100.0 * hist[i]->Integral(integralBegin, integralEnd) / hist[i]->GetEntries()) : -1.0);
1697     Printf("");
1698   }
1699 }
1700
1701 TH2* GetCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
1702 {
1703   // returns the correction factor with pt integrated out
1704
1705   loadlibs();
1706
1707   TFile::Open(fileName);
1708
1709   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
1710   if (!dNdEtaCorrection->LoadHistograms())
1711     return;
1712
1713   //  hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
1714
1715   gener = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1716   measu = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1717
1718   gener->GetZaxis()->SetRange(gener->GetZaxis()->FindBin(ptmin), gener->GetNbinsZ()+1);
1719   TH2D *gener_xy = gener->Project3D("yx");
1720
1721   measu->GetZaxis()->SetRange(measu->GetZaxis()->FindBin(ptmin), measu->GetNbinsZ()+1);
1722   TH2D *measu_xy = measu->Project3D("yx");
1723
1724   cout << measu->GetZaxis()->FindBin(ptmin) << " " << measu->GetNbinsZ()+1 << endl;
1725
1726   TCanvas *canp = new TCanvas("canp","canp",600,1000);
1727   canp->Divide(1,2,0.0001,0.0001);
1728   canp->cd(1);
1729   gener_xy->Draw("COLZ");
1730   canp->cd(2);
1731   measu_xy->Draw("COLZ");
1732
1733
1734   TCanvas *canpr = new TCanvas("canpr","canpr",700,500);
1735   canpr->cd();
1736   TH2D *proj = new TH2D(*gener_xy);
1737   proj->Divide(measu_xy);
1738
1739 //   proj = hist->Project3D("yx");
1740   proj->Draw("COLZ");
1741
1742   return proj;
1743 }
1744
1745 void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
1746 {
1747   TH2* proj = GetCorrection(fileName, dirName, ptmin);
1748
1749   const Float_t limit = 5;
1750
1751   TString array = "{";
1752   TString arrayEnd = "}";
1753
1754   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
1755   {
1756     Int_t begin = -1;
1757     Int_t end = -1;
1758     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
1759     {
1760       if (begin == -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
1761         begin = x;
1762       if (begin != -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
1763         end = x;
1764     }
1765     Printf("Limits for y = %d are %d to %d", y, begin, end);
1766
1767     if (y > 1)
1768       array += ", ";
1769     array += Form("%d", begin);
1770
1771     if (y > 1)
1772       arrayEnd.Prepend(", ");
1773     arrayEnd.Prepend(Form("%d", (end == -1) ? -1 : proj->GetNbinsX() + 1 - end));
1774   }
1775   array += "}";
1776   arrayEnd.Prepend("{");
1777
1778   Printf("Begin array:");
1779   Printf("%s", array.Data());
1780
1781   Printf("End array (mirrored) (should be the same):");
1782   Printf("%s", arrayEnd.Data());
1783 }
1784
1785 void AverageMultiplicity(const char* fileName = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1786 {
1787   loadlibs();
1788
1789   TFile::Open(fileName);
1790
1791   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1792   dNdEtaCorrection->LoadHistograms();
1793   TH2* events = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram();
1794   TH3* tracks = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram();
1795
1796   Float_t nEvents = events->Integral(events->GetXaxis()->FindBin(-1), events->GetXaxis()->FindBin(1), 0, events->GetNbinsY()+1);
1797   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);
1798
1799   Printf("%f %f --> %f", nEvents, nTracks, nTracks / nEvents);
1800 }
1801
1802 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")
1803 {
1804   loadlibs();
1805
1806   TFile::Open(rawFile);
1807   dNdEtaAnalysis* raw = new dNdEtaAnalysis("dndeta", "dndeta");
1808   raw->LoadHistograms("fdNdEtaAnalysisESD");
1809   raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->GetXaxis()->SetRangeUser(-vertexRange, vertexRange);
1810   tracks = raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Project3D("y");
1811   events = raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("events", 0, raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->GetNbinsY() + 1);
1812   Float_t nEvents = events->Integral(events->FindBin(-vertexRange), events->FindBin(vertexRange));
1813   tracks->Scale(1.0 / nEvents / tracks->GetBinWidth(1));
1814
1815   TFile::Open(mcFile);
1816   dNdEtaAnalysis* mc = new dNdEtaAnalysis("dndeta", "dndeta");
1817   mc->LoadHistograms("dndetaTrVtx");
1818   mcH = mc->GetdNdEtaPtCutOffCorrectedHistogram(0);
1819
1820   new TCanvas;
1821   mcH->SetLineColor(2);
1822   mcH->DrawCopy();
1823   tracks->DrawCopy("SAME");
1824
1825   new TCanvas;
1826   mcH->GetYaxis()->SetRangeUser(0, 5);
1827   mcH->Divide(tracks);
1828   mcH->DrawCopy();
1829   mcH->Fit("pol0", "", "", -etaRange, etaRange);
1830 }
1831
1832 void TrackCuts_Comparison(char* histName, Int_t plotWhich = 0, const char* fileName = "correction_map.root")
1833 {
1834   // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
1835   //    --> manually disable it in the run.C
1836   //
1837   // plotWhich: 0 = only before
1838   //            1 = both
1839   //            2 = only after
1840
1841   file = TFile::Open(fileName);
1842
1843   Int_t count = 0;
1844   Int_t colors[] = { 1, 2, 3, 4, 5, 6 };
1845
1846   TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
1847   TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
1848   TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
1849
1850   TCanvas* c1 = new TCanvas("c1", "c1", 800, 1200);
1851   c1->Divide(1, 2);
1852   //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
1853   //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
1854
1855   const char* folders2[] = { "before_cuts", "after_cuts" };
1856   Bool_t first = kTRUE;
1857   for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
1858   {
1859     const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
1860     const char* names[] =    { "all", "primaries", "secondaries" };
1861     TH1* base = 0;
1862     TH1* prim = 0;
1863     TH1* sec = 0;
1864     for (Int_t i = 0; i < 3; i++)
1865     {
1866       TString folder;
1867       folder.Form("%s/%s/%s", folders1[i], folders2[j], histName);
1868       TH1* hist = (TH1*) file->Get(folder);
1869       legend->AddEntry(hist, Form("%s %s", names[i], folders2[j]));
1870
1871       c1->cd(1);
1872       hist->SetLineColor(colors[count]);
1873       hist->DrawCopy((count == 0) ? "" : "SAME");
1874
1875       switch (i)
1876       {
1877         case 0: base = hist; break;
1878         case 1: prim = hist; break;
1879         case 2: sec = hist; break;
1880       }
1881
1882       count++;
1883     }
1884
1885     TH1* eff    = (TH1*) prim->Clone("eff"); eff->Reset();
1886     TH1* purity = (TH1*) prim->Clone("purity"); purity->Reset();
1887
1888     for (Int_t bin = 1; bin <= prim->GetNbinsX(); bin++)
1889     {
1890       eff->SetBinContent(bin, prim->Integral(1, bin) / prim->Integral(1, prim->GetNbinsX() + 1));
1891       if (prim->Integral(1, bin) + sec->Integral(1, bin) > 0)
1892         purity->SetBinContent(bin, sec->Integral(1, bin) / (prim->Integral(1, bin) + sec->Integral(1, bin)));
1893     }
1894
1895     eff->GetYaxis()->SetRangeUser(0, 1);
1896     eff->SetLineColor(colors[0+j*2]);
1897     eff->SetStats(kFALSE);
1898     purity->SetLineColor(colors[1+j*2]);
1899
1900     legend3->AddEntry(eff, Form("%s: efficiency", folders2[j]));
1901     legend3->AddEntry(purity, Form("%s: contamination", folders2[j]));
1902
1903     c1->cd(2);
1904     eff->DrawCopy((first) ? "" : "SAME");
1905     first = kFALSE;
1906     purity->DrawCopy("SAME");
1907   }
1908
1909   c1->cd(1)->SetLogy();
1910   c1->cd(1)->SetGridx();
1911   c1->cd(1)->SetGridy();
1912   legend->Draw();
1913
1914   //c2->cd();
1915  // c2->SetGridx();
1916  // c2->SetGridy();
1917   //legend2->Draw();
1918
1919   c1->cd(2)->SetGridx();
1920   c1->cd(2)->SetGridy();
1921   legend3->Draw();
1922
1923   c1->SaveAs(Form("%s.png", histName));
1924 }
1925
1926 void TrackCuts_DCA()
1927 {
1928   file = TFile::Open("correction_map.root");
1929   hist = (TH2*) file->Get("esd_track_cuts/before_cuts/dXYvsDZ");
1930
1931   TCanvas* c1 = new TCanvas("c1", "c1", 600, 600);
1932   c1->SetLogz();
1933   c1->SetRightMargin(0.12);
1934   c1->SetBottomMargin(0.12);
1935
1936   hist->SetStats(kFALSE);
1937   hist->Draw("COLZ");
1938
1939   ellipse = new TEllipse(0, 0, 4);
1940   ellipse->SetLineWidth(2);
1941   ellipse->SetLineStyle(2);
1942   ellipse->SetFillStyle(0);
1943   ellipse->Draw();
1944
1945   c1->SaveAs("trackcuts_dca_2d.eps");
1946 }
1947
1948 void FindNSigma(TH2* hist, Int_t nSigma = 1)
1949 {
1950   TH1* proj = hist->ProjectionY();
1951   proj->Reset();
1952
1953   for (Int_t bin=1; bin<=proj->GetNbinsX(); bin++)
1954   {
1955     if (hist->Integral(1, hist->GetNbinsX(), bin, bin) == 0)
1956       continue;
1957
1958     Int_t limit = -1;
1959     for (limit = 1; limit<=hist->GetNbinsX(); limit++)
1960   }
1961 }
1962
1963 void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
1964 {
1965   TH2* proj = GetCorrection(fileName, dirName, ptmin);
1966
1967   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
1968     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
1969       if (proj->GetBinContent(x, y) > 5 || proj->GetBinContent(x, y) == 0)
1970       {
1971         proj->SetBinContent(x, y, 0);
1972       }
1973       else
1974         proj->SetBinContent(x, y, 1);
1975
1976
1977   input->Multiply(proj);
1978 }