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