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