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