merging post-CVS changes into SVN:
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / drawPlots.C
1 /* $Id$ */
2
3 #if !defined(__CINT__) || defined(__MAKECINT__)
4
5 #include "AliPWG0Helper.h"
6 #include "dNdEtaAnalysis.h"
7 #include "AlidNdEtaCorrection.h"
8
9 #include <TCanvas.h>
10 #include <TFile.h>
11 #include <TH1.h>
12 #include <TH2F.h>
13 #include <TH3F.h>
14 #include <TLine.h>
15 #include <TSystem.h>
16
17 #endif
18
19 Int_t gMax = 5;
20
21 extern TSystem* gSystem;
22
23 void loadlibs()
24 {
25   gSystem->Load("libANALYSIS");
26   gSystem->Load("libPWG0base");
27 }
28
29 void SetRanges(TAxis* axis)
30 {
31   if (strcmp(axis->GetTitle(), "#eta") == 0)
32     axis->SetRangeUser(-1.7999, 1.7999);
33   if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
34     axis->SetRangeUser(0, 9.9999);
35   if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
36     axis->SetRangeUser(-15, 14.9999);
37   if (strcmp(axis->GetTitle(), "Ntracks") == 0)
38     axis->SetRangeUser(0, 99.9999);
39 }
40
41 void SetRanges(TH1* hist)
42 {
43   SetRanges(hist->GetXaxis());
44   SetRanges(hist->GetYaxis());
45   SetRanges(hist->GetZaxis());
46 }
47
48
49 void Prepare3DPlot(TH3* hist)
50 {
51   hist->GetXaxis()->SetTitleOffset(1.5);
52   hist->GetYaxis()->SetTitleOffset(1.5);
53   hist->GetZaxis()->SetTitleOffset(1.5);
54
55   hist->SetStats(kFALSE);
56 }
57
58 void Prepare2DPlot(TH2* hist)
59 {
60   hist->SetStats(kFALSE);
61   hist->GetYaxis()->SetTitleOffset(1.4);
62
63   hist->SetMinimum(0);
64   hist->SetMaximum(gMax);
65
66   SetRanges(hist);
67 }
68
69 void Prepare1DPlot(TH1* hist)
70 {
71   hist->SetLineWidth(2);
72   hist->SetStats(kFALSE);
73
74   hist->GetXaxis()->SetLabelOffset(0.02);
75   hist->GetXaxis()->SetTitleOffset(1.3);
76   hist->GetYaxis()->SetTitleOffset(1.3);
77
78   SetRanges(hist);
79 }
80
81 void InitPad()
82 {
83   gPad->Range(0, 0, 1, 1);
84   gPad->SetLeftMargin(0.15);
85   //gPad->SetRightMargin(0.05);
86   //gPad->SetTopMargin(0.13);
87   gPad->SetBottomMargin(0.12);
88
89   gPad->SetGridx();
90   gPad->SetGridy();
91 }
92
93 void InitPadCOLZ()
94 {
95   gPad->Range(0, 0, 1, 1);
96   gPad->SetRightMargin(0.15);
97   gPad->SetLeftMargin(0.12);
98
99   gPad->SetGridx();
100   gPad->SetGridy();
101 }
102
103 // --- end of helpers --- begin functions ---
104
105 void DrawOverview(const char* fileName, const char* dirName)
106 {
107   loadlibs();
108   if (!TFile::Open(fileName))
109     return;
110
111   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
112   if (!dNdEtaCorrection->LoadHistograms())
113     return;
114
115   dNdEtaCorrection->Finish();
116
117   dNdEtaCorrection->DrawOverview();
118 }
119
120 void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
121 {
122   gSystem->Load("libPWG0base");
123
124   TFile::Open(esdFile);
125   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
126   fdNdEtaAnalysisESD->LoadHistograms();
127
128   TFile::Open(mcFile);
129   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
130   fdNdEtaAnalysisMC->LoadHistograms();
131   //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
132
133   for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
134     fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
135
136   fdNdEtaAnalysisESD->DrawHistograms();
137 }
138
139 void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
140 {
141   gSystem->Load("libPWG0base");
142
143   const char* ESDfolder = 0;
144
145   if (plot == 0) // all
146     ESDfolder = "dndeta";
147   else if (plot == 1) // mb
148     ESDfolder = "dndeta_mb";
149   else if (plot == 2) // mb vtx
150     ESDfolder = "dndeta_mbvtx";
151
152   TFile::Open("analysis_esd.root");
153   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
154   fdNdEtaAnalysisESD->LoadHistograms();
155
156   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
157   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
158
159   TH2F* hist = 0;
160
161   if (plot == 0) // all
162     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
163   else if (plot == 1) // mb
164     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
165   else if (plot == 2) // mb vtx
166     hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
167
168   TH1* proj = hist->ProjectionX();
169
170   TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
171   for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
172   {
173     Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
174     if (value != 0)
175     {
176       printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
177       vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
178     }
179   }
180
181   new TCanvas;
182   vertex->DrawCopy();
183 }
184
185 void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
186 {
187   gSystem->Load("libPWG0base");
188
189   TFile::Open("analysis_esd.root");
190   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
191   fdNdEtaAnalysisESD->LoadHistograms();
192
193   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
194   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
195
196   //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
197   //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
198
199   TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
200   TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
201
202   TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
203
204   new TCanvas;
205   histESD->Draw();
206
207   new TCanvas;
208   histCorr->Draw();
209
210   for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
211     for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
212       for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
213       {
214         Float_t value1 = histESD->GetBinContent(x, y, z);
215         Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
216
217         if (value2 > 0 && value1 > 0)
218         {
219           printf("%f %f %f\n", value1, value2, value1 / value2);
220           diff->Fill(value1 / value2);
221         }
222       }
223
224   new TCanvas;
225   diff->Draw();
226 }
227
228 Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
229 {
230   Double_t avgMC = 0;
231   Double_t avgESD = 0;
232   for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
233   {
234     avgMC += histMC->GetBinContent(bin);
235     avgESD += histESD->GetBinContent(bin);
236   }
237   Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
238
239   avgMC /= nBins;
240   avgESD /= nBins;
241
242   // deviation when integrate in |eta| < 0.8 between mc and esd
243   Double_t diffFullRange = (avgMC - avgESD) / avgMC;
244
245   Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
246
247   return diffFullRange;
248 }
249
250 void dNdEta(Bool_t onlyESD = kFALSE)
251 {
252   TFile* file = TFile::Open("analysis_esd.root");
253   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
254   TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
255   TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
256   TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
257   TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
258   TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
259   TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
260
261   Prepare1DPlot(histESD);
262   Prepare1DPlot(histESDMB);
263   Prepare1DPlot(histESDMBVtx);
264
265   Prepare1DPlot(histESDNoPt);
266   Prepare1DPlot(histESDMBNoPt);
267   Prepare1DPlot(histESDMBVtxNoPt);
268   Prepare1DPlot(histESDMBTracksNoPt);
269
270   histESD->SetLineWidth(0);
271   histESDMB->SetLineWidth(0);
272   histESDMBVtx->SetLineWidth(0);
273
274   histESDNoPt->SetLineWidth(0);
275   histESDMBNoPt->SetLineWidth(0);
276   histESDMBVtxNoPt->SetLineWidth(0);
277
278   histESD->SetMarkerColor(1);
279   histESDMB->SetMarkerColor(2);
280   histESDMBVtx->SetMarkerColor(3);
281
282   histESDNoPt->SetMarkerColor(1);
283   histESDMBNoPt->SetMarkerColor(2);
284   histESDMBVtxNoPt->SetMarkerColor(3);
285   histESDMBTracksNoPt->SetMarkerColor(4);
286
287   histESD->SetMarkerStyle(20);
288   histESDMB->SetMarkerStyle(21);
289   histESDMBVtx->SetMarkerStyle(22);
290
291   histESDNoPt->SetMarkerStyle(20);
292   histESDMBNoPt->SetMarkerStyle(21);
293   histESDMBVtxNoPt->SetMarkerStyle(22);
294   histESDMBTracksNoPt->SetMarkerStyle(23);
295
296   TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 1000, 0, histESDMBVtx->GetMaximum() * 1.1);
297   Prepare1DPlot(dummy);
298   dummy->SetStats(kFALSE);
299   dummy->SetXTitle("#eta");
300   dummy->SetYTitle("dN_{ch}/d#eta");
301   dummy->GetYaxis()->SetTitleOffset(1);
302
303   Float_t etaLimit = 0.7999;
304
305   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
306   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
307   histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
308
309   histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
310   histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
311   histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
312   histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
313
314   TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
315
316   dummy->DrawCopy();
317   histESDMBVtx->Draw("SAME");
318   histESDMB->Draw("SAME");
319   histESD->Draw("SAME");
320
321   canvas->SaveAs("dNdEta1.gif");
322   canvas->SaveAs("dNdEta1.eps");
323
324   if (onlyESD)
325     return;
326
327   loadlibs();
328
329   TFile* file2 = TFile::Open("analysis_mc.root");
330   TH1* histMC = (TH1*) file2->Get("dndeta/dNdEta_corrected")->Clone("cloned");
331   TH1* histMCTr = (TH1*) file2->Get("dndetaTr/dNdEta_corrected")->Clone("cloned2");
332   TH1* histMCTrVtx = (TH1*) file2->Get("dndetaTrVtx/dNdEta_corrected")->Clone("cloned3");
333
334   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
335   fdNdEtaAnalysis->LoadHistograms();
336   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
337   TH1* histMCPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
338
339   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
340   fdNdEtaAnalysis->LoadHistograms();
341   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
342   TH1* histMCTrPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
343
344   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
345   fdNdEtaAnalysis->LoadHistograms();
346   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
347   TH1* histMCTrVtxPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
348
349   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
350   fdNdEtaAnalysis->LoadHistograms();
351   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
352   TH1* histMCTracksPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
353
354   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
355
356   Prepare1DPlot(histMC);
357   Prepare1DPlot(histMCTr);
358   Prepare1DPlot(histMCTrVtx);
359
360   Prepare1DPlot(histMCPtCut);
361   Prepare1DPlot(histMCTrPtCut);
362   Prepare1DPlot(histMCTrVtxPtCut);
363   if (histMCTracksPtCut)
364     Prepare1DPlot(histMCTracksPtCut);
365
366   histMC->SetLineColor(1);
367   histMCTr->SetLineColor(2);
368   histMCTrVtx->SetLineColor(3);
369
370   histMCPtCut->SetLineColor(1);
371   histMCTrPtCut->SetLineColor(2);
372   histMCTrVtxPtCut->SetLineColor(3);
373   if (histMCTracksPtCut)
374     histMCTracksPtCut->SetLineColor(4);
375
376   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
377   Prepare1DPlot(dummy2);
378   dummy2->GetYaxis()->SetRangeUser(0, histESDMBVtx->GetMaximum() * 1.1);
379
380   dummy2->DrawCopy();
381   histMC->Draw("SAME");
382   histMCTr->Draw("SAME");
383   histMCTrVtx->Draw("SAME");
384   histESD->Draw("SAME");
385   histESDMB->Draw("SAME");
386   histESDMBVtx->Draw("SAME");
387   histESDNoPt->Draw("SAME");
388   histESDMBNoPt->Draw("SAME");
389   histESDMBVtxNoPt->Draw("SAME");
390   histESDMBTracksNoPt->Draw("SAME");
391   histMCPtCut->Draw("SAME");
392   histMCTrPtCut->Draw("SAME");
393   histMCTrVtxPtCut->Draw("SAME");
394   if (histMCTracksPtCut)
395     histMCTracksPtCut->Draw("SAME");
396
397   canvas2->SaveAs("dNdEta2.gif");
398   canvas2->SaveAs("dNdEta2.eps");
399
400   TH1* ratio = (TH1*) histMC->Clone("ratio");
401   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
402
403   ratio->Divide(histESD);
404   ratioNoPt->Divide(histESDNoPt);
405
406   ratio->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
407
408   ratio->SetLineColor(1);
409   ratioNoPt->SetLineColor(2);
410
411   Double_t average = 0;       // average deviation from 1 in ratio (depends on the number of bins if statistical)
412   for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
413     average += TMath::Abs(ratio->GetBinContent(bin) - 1);
414   Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
415   average /= nBins;
416   Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
417
418   PrintIntegratedDeviation(histMC, histESD, "all events");
419   PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
420   PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
421   PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
422   PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
423   PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
424
425   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 600);
426   canvas3->Range(0, 0, 1, 1);
427   //canvas3->Divide(1, 2, 0, 0);
428
429   //canvas3->cd(1);
430   TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
431   pad1->Draw();
432
433   TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
434   pad2->Draw();
435
436   pad1->SetRightMargin(0.05);
437   pad2->SetRightMargin(0.05);
438
439   // no border between them
440   pad1->SetBottomMargin(0);
441   pad2->SetTopMargin(0);
442
443   pad1->cd();
444
445   TLegend* legend = new TLegend(0.4, 0.05, 0.65, 0.3);
446   legend->SetFillColor(0);
447   legend->AddEntry(histESDMBVtx, "triggered, vertex");
448   legend->AddEntry(histESDMB, "triggered");
449   legend->AddEntry(histESD, "all events");
450   legend->AddEntry(histMC, "MC prediction");
451
452   dummy->GetXaxis()->SetLabelSize(0.06);
453   dummy->GetYaxis()->SetLabelSize(0.06);
454   dummy->GetXaxis()->SetTitleSize(0.06);
455   dummy->GetYaxis()->SetTitleSize(0.06);
456   dummy->GetYaxis()->SetTitleOffset(0.7);
457   dummy->DrawCopy();
458   histESDMBVtx->Draw("SAME");
459   histESDMB->Draw("SAME");
460   histESD->Draw("SAME");
461   histMC->Draw("SAME");
462
463   legend->Draw();
464
465   pad2->cd();
466   pad2->SetBottomMargin(0.15);
467
468   Float_t min = TMath::Min(0.961, ratio->GetMinimum() * 0.95);
469   Float_t max = TMath::Max(1.049, ratio->GetMaximum() * 1.05);
470
471   TH1F dummy3("dummy3", ";#eta;Ratio: MC / ESD", 1, -1.5, 1.5);
472   dummy3.SetStats(kFALSE);
473   dummy3.SetBinContent(1, 1);
474   dummy3.GetYaxis()->SetRangeUser(min, max);
475   dummy3.SetLineWidth(2);
476   dummy3.GetXaxis()->SetLabelSize(0.06);
477   dummy3.GetYaxis()->SetLabelSize(0.06);
478   dummy3.GetXaxis()->SetTitleSize(0.06);
479   dummy3.GetYaxis()->SetTitleSize(0.06);
480   dummy3.GetYaxis()->SetTitleOffset(0.7);
481   dummy3.DrawCopy();
482
483   ratio->Draw("SAME");
484
485   //pad2->Draw();
486
487   canvas3->Modified();
488
489   canvas3->SaveAs("dNdEta.gif");
490   canvas3->SaveAs("dNdEta.eps");
491
492   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
493
494   ratio->Draw();
495   ratioNoPt->Draw("SAME");
496
497   TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
498   legend->SetFillColor(0);
499   legend->AddEntry(ratio, "mc/esd");
500   legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
501   legend->Draw();
502 }
503
504 void ptSpectrum()
505 {
506   TFile* file = TFile::Open("analysis_esd.root");
507   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
508
509   TFile* file2 = TFile::Open("analysis_mc.root");
510   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
511
512   TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
513   InitPad();
514   gPad->SetLogy();
515
516   Prepare1DPlot(histMC);
517   Prepare1DPlot(histESD);
518
519   histESD->SetTitle("");
520   histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
521   histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
522
523   histMC->SetLineColor(kBlue);
524   histESD->SetLineColor(kRed);
525
526   histESD->GetYaxis()->SetTitleOffset(1.5);
527   histESD->GetXaxis()->SetRangeUser(0, 4.9999);
528
529   histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
530
531   histESD->Draw();
532   histMC->Draw("SAME");
533
534   canvas->SaveAs("ptSpectrum.gif");
535   canvas->SaveAs("ptSpectrum.eps");
536 }
537
538 void ptCutoff()
539 {
540   gSystem->Load("libPWG0base");
541
542   TFile::Open("correction_map.root");
543   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
544   dNdEtaCorrection->LoadHistograms();
545
546   dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, -100, kTRUE);
547
548   TH1* hist = dynamic_cast<TH1*> (gROOT->FindObject("generated_pt")->Clone("ptcutoff"));
549
550   hist->GetXaxis()->SetRangeUser(0, 0.9999);
551   hist->SetMinimum(0);
552
553   hist->SetTitle("Generated Particles");
554   Prepare1DPlot(hist);
555
556   TCanvas* canvas = new TCanvas("ptCutoff", "ptCutoff", 700, 500);
557   hist->DrawCopy();
558
559   TLine* line = new TLine(0.3, 0 - hist->GetMaximum() * 0, 0.3, hist->GetMaximum() * 1.1);
560   line->SetLineWidth(3);
561   line->SetLineColor(kRed);
562   line->Draw();
563
564   canvas->SaveAs("ptCutoff.gif");
565   canvas->SaveAs("ptCutoff.eps");
566
567   TH1F* factor = new TH1F("factor", ";#eta;correction factor", 20, -1, 1.000001);
568   factor->SetLineWidth(2);
569   for (Float_t eta = -0.95; eta<1; eta += 0.1)
570     factor->Fill(eta, 1.0 / dNdEtaCorrection->GetMeasuredFraction(AlidNdEtaCorrection::kINEL, 0.3, eta, kFALSE));
571
572   TCanvas* canvas = new TCanvas("ptCutoff_factor", "ptCutoff_factor", 700, 500);
573   InitPad();
574
575   Prepare1DPlot(factor);
576   factor->GetYaxis()->SetRangeUser(1, 2);
577   factor->GetYaxis()->SetTitleOffset(1);
578   factor->Draw();
579
580   canvas->SaveAs("ptCutoff_factor.eps");
581 }
582
583 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
584 {
585   gSystem->Load("libPWG0base");
586
587   TFile::Open(fileName);
588   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
589   dNdEtaCorrection->LoadHistograms();
590
591   TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
592   TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
593
594   Prepare2DPlot(corrTrigger);
595   corrTrigger->SetTitle("b) Trigger bias correction");
596
597   Prepare2DPlot(corrVtx);
598   corrVtx->SetTitle("a) Vertex reconstruction correction");
599
600   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
601   corrVtx->GetYaxis()->SetTitle("Multiplicity");
602
603   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
604   canvas->Divide(2, 1);
605
606   canvas->cd(1);
607   InitPadCOLZ();
608   corrVtx->DrawCopy("COLZ");
609
610   canvas->cd(2);
611   InitPadCOLZ();
612   corrTrigger->DrawCopy("COLZ");
613
614   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
615   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
616
617   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
618   canvas->Divide(2, 1);
619
620   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
621   corrVtx->GetYaxis()->SetRangeUser(0, 5);
622
623   canvas->cd(1);
624   InitPadCOLZ();
625   corrVtx->DrawCopy("COLZ");
626
627   canvas->cd(2);
628   InitPadCOLZ();
629   corrTrigger->DrawCopy("COLZ");
630
631   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
632   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
633 }
634
635 void TriggerBias(const char* fileName = "correction_map.root")
636 {
637   TFile* file = TFile::Open(fileName);
638
639   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
640
641   Prepare2DPlot(corr);
642   corr->SetTitle("Trigger bias correction");
643
644   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
645   InitPadCOLZ();
646   corr->DrawCopy("COLZ");
647
648   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
649   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
650
651   corr->GetYaxis()->SetRangeUser(0, 5);
652
653   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
654   InitPadCOLZ();
655   corr->DrawCopy("COLZ");
656
657   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
658   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
659 }
660
661 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
662 {
663   gSystem->Load("libPWG0base");
664
665   TFile* file = TFile::Open(fileName);
666   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
667   dNdEtaCorrection->LoadHistograms();
668
669   TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
670   TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
671
672   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
673   canvas->Divide(2, 1);
674
675   canvas->cd(1);
676   InitPad();
677
678   Prepare1DPlot(hist);
679   hist->SetTitle("");
680   hist->GetYaxis()->SetTitle("correction factor");
681   hist->GetYaxis()->SetRangeUser(1, 1.5);
682   hist->GetYaxis()->SetTitleOffset(1.6);
683   hist->Draw();
684
685   canvas->cd(2);
686   InitPad();
687
688   Prepare1DPlot(hist2);
689   hist2->SetTitle("");
690   hist2->GetYaxis()->SetTitle("correction factor");
691   hist2->GetXaxis()->SetRangeUser(0, 5);
692   hist2->GetYaxis()->SetTitleOffset(1.6);
693   hist2->GetXaxis()->SetTitle("multiplicity");
694   hist2->Draw();
695
696   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
697   pave->SetFillColor(0);
698   pave->AddText("|z| < 10 cm");
699   pave->Draw();
700
701   canvas->SaveAs("TriggerBias1D.eps");
702 }
703
704 void VtxRecon()
705 {
706   TFile* file = TFile::Open("correction_map.root");
707
708   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
709
710   Prepare2DPlot(corr);
711   corr->SetTitle("Vertex reconstruction correction");
712
713   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
714   InitPadCOLZ();
715   corr->DrawCopy("COLZ");
716
717   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
718   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
719
720   corr->GetYaxis()->SetRangeUser(0, 5);
721
722   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
723   InitPadCOLZ();
724   corr->DrawCopy("COLZ");
725
726   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
727   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
728 }
729
730 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
731 {
732   gSystem->Load("libPWG0base");
733
734   TFile* file = TFile::Open(fileName);
735   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
736   dNdEtaCorrection->LoadHistograms();
737
738   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
739   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
740
741   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
742   canvas->Divide(2, 1);
743
744   canvas->cd(1);
745   InitPad();
746
747   Prepare1DPlot(hist);
748   hist->SetTitle("");
749   hist->GetYaxis()->SetTitle("correction factor");
750   hist->GetYaxis()->SetRangeUser(1, 1.8);
751   hist->GetYaxis()->SetTitleOffset(1.6);
752   hist->DrawCopy();
753
754   canvas->cd(2);
755   InitPad();
756
757   Prepare1DPlot(hist2);
758   hist2->SetTitle("");
759   hist2->GetYaxis()->SetTitle("correction factor");
760   hist2->GetXaxis()->SetRangeUser(0, 20);
761   hist2->GetYaxis()->SetTitleOffset(1.6);
762   hist2->GetXaxis()->SetTitle("multiplicity");
763   hist2->Draw();
764
765   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
766   pave->SetFillColor(0);
767   pave->AddText("|z| < 10 cm");
768   pave->Draw();
769
770   canvas->SaveAs("VtxRecon1D.eps");
771
772   Correction1DCreatePlots(fileName, folderName, 9.9, 2);
773
774   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
775   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
776
777   Prepare1DPlot(corrX);
778   Prepare1DPlot(corrZ);
779
780   corrX->GetYaxis()->SetTitleOffset(1.5);
781   corrZ->GetYaxis()->SetTitleOffset(1.5);
782
783   corrX->SetTitle("a) z projection");
784   corrZ->SetTitle("b) p_{T} projection");
785
786   corrX->GetYaxis()->SetTitle("correction factor");
787   corrZ->GetYaxis()->SetTitle("correction factor");
788
789   corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
790
791   TString canvasName;
792   canvasName.Form("VtxRecon1D_Track");
793   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
794   canvas->Divide(2, 1);
795
796   canvas->cd(1);
797   InitPad();
798   corrX->DrawCopy();
799
800   canvas->cd(2);
801   InitPad();
802   gPad->SetLogx();
803   corrZ->Draw();
804
805   canvas->SaveAs("VtxRecon1D_Track.eps");
806   canvas->SaveAs("VtxRecon1D_Track.gif");
807 }
808
809 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
810 {
811   gSystem->Load("libPWG0base");
812
813   TFile::Open(fileName);
814   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
815   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
816
817   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
818   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
819
820   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
821   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
822
823   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
824   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
825   gene->GetXaxis()->SetRangeUser(-10, 10);
826   meas->GetXaxis()->SetRangeUser(-10, 10);
827
828   Float_t eff1 = gene->Integral() / meas->Integral();
829   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
830
831   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
832
833   gene->GetZaxis()->SetRangeUser(0.3, 10);
834   meas->GetZaxis()->SetRangeUser(0.3, 10);
835
836   Float_t eff2 = gene->Integral() / meas->Integral();
837   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
838
839   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
840
841   gene->GetZaxis()->SetRangeUser(0.3, 1);
842   meas->GetZaxis()->SetRangeUser(0.3, 1);
843
844   Float_t eff3 = gene->Integral() / meas->Integral();
845   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
846
847   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
848 }
849
850 void Correction1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType = 0)
851 {
852   TFile::Open(fileName);
853   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
854   dNdEtaCorrection->LoadHistograms();
855
856   TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
857   TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetMeasuredHistogram();
858
859   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
860   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
861   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
862   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
863   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kTRUE);
864   gene->GetYaxis()->SetRange(0, 0);
865   meas->GetYaxis()->SetRange(0, 0);
866
867   gene->GetXaxis()->SetRangeUser(-10, 10);
868   meas->GetXaxis()->SetRangeUser(-10, 10);
869   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kTRUE);
870   gene->GetZaxis()->SetRange(0, 0);
871   meas->GetZaxis()->SetRange(0, 0);
872
873   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
874   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
875   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kTRUE);
876 }
877
878 void Correction1D(Int_t correctionType = 0, const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
879 {
880   gSystem->Load("libPWG0base");
881
882   Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType);
883
884   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
885   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
886   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
887
888   Prepare1DPlot(corrX);
889   Prepare1DPlot(corrY);
890   Prepare1DPlot(corrZ);
891
892   corrX->SetTitle("a) z projection");
893   corrY->SetTitle("b) #eta projection");
894   corrZ->SetTitle("c) p_{T} projection");
895
896   corrX->GetYaxis()->SetTitle("correction factor");
897   corrY->GetYaxis()->SetTitle("correction factor");
898   corrZ->GetYaxis()->SetTitle("correction factor");
899
900   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
901
902   TString canvasName;
903   canvasName.Form("Correction1D_%s", folder);
904   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
905   canvas->Divide(3, 1);
906
907   canvas->cd(1);
908   InitPad();
909   corrX->DrawCopy();
910
911   canvas->cd(2);
912   InitPad();
913   corrY->Draw();
914
915   canvas->cd(3);
916   InitPad();
917   corrZ->Draw();
918
919   canvas->SaveAs(Form("Correction1D_%d_%s_%f.gif", correctionType, fileName, upperPtLimit));
920   canvas->SaveAs(Form("Correction1D_%d_%s_%f.eps", correctionType, fileName, upperPtLimit));
921 }
922
923 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
924 {
925   gSystem->Load("libPWG0base");
926
927   Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
928
929   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
930   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
931   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
932
933   Prepare1DPlot(corrX);
934   Prepare1DPlot(corrY);
935   Prepare1DPlot(corrZ);
936
937   corrX->SetTitle("a) z projection");
938   corrY->SetTitle("a) #eta projection");
939   corrZ->SetTitle("b) p_{T} projection");
940
941   corrY->GetYaxis()->SetTitle("correction factor");
942   corrZ->GetYaxis()->SetTitle("correction factor");
943
944   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
945
946   TString canvasName;
947   canvasName.Form("Track2Particle1D_%s", folder);
948   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
949   canvas->Divide(3, 1);
950
951   canvas->cd(1);
952   InitPad();
953   corrX->DrawCopy();
954
955   canvas->cd(2);
956   InitPad();
957   corrY->Draw();
958
959   canvas->cd(3);
960   InitPad();
961   corrZ->Draw();
962
963   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
964   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
965
966   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
967
968   canvasName.Form("Track2Particle1D_%s_etapt", folder);
969   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
970   canvas->Divide(2, 1);
971
972   canvas->cd(1);
973   InitPad();
974   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
975   corrY->GetYaxis()->SetRangeUser(1, 1.5);
976   corrY->GetYaxis()->SetTitleOffset(1.5);
977   corrY->DrawCopy();
978   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
979   pave->AddText("|z| < 10 cm");
980   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
981   pave->Draw();
982
983   canvas->cd(2);
984   InitPad();
985   gPad->SetLogx();
986   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
987   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
988   corrZ->GetYaxis()->SetTitleOffset(1.5);
989   corrZ->DrawCopy();
990   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
991   pave->AddText("|z| < 10 cm");
992   pave->AddText("|#eta| < 0.8");
993   pave->Draw();
994
995   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
996   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
997 }
998
999 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1000 {
1001   gSystem->Load("libPWG0base");
1002
1003   // particle type
1004   for (Int_t particle=0; particle<4; ++particle)
1005   {
1006     TString dirName;
1007     dirName.Form("correction_%d", particle);
1008     Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
1009
1010     TString tmpx, tmpy, tmpz;
1011     tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1012     tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1013     tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
1014
1015     TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1016     TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1017     TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
1018
1019     Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
1020
1021     TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1022     TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1023     TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
1024
1025     posX->Divide(negX);
1026     posY->Divide(negY);
1027     posZ->Divide(negZ);
1028
1029     Prepare1DPlot(posX);
1030     Prepare1DPlot(posY);
1031     Prepare1DPlot(posZ);
1032
1033     Float_t min = 0.8;
1034     Float_t max = 1.2;
1035
1036     posX->SetMinimum(min);
1037     posX->SetMaximum(max);
1038     posY->SetMinimum(min);
1039     posY->SetMaximum(max);
1040     posZ->SetMinimum(min);
1041     posZ->SetMaximum(max);
1042
1043     posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1044
1045     posX->GetYaxis()->SetTitleOffset(1.7);
1046     posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1047     posY->GetYaxis()->SetTitleOffset(1.7);
1048     posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1049     posZ->GetYaxis()->SetTitleOffset(1.7);
1050     posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
1051
1052     posZ->GetXaxis()->SetRangeUser(0, 1);
1053
1054     TString canvasName;
1055     canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
1056
1057     TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1058     canvas->Divide(3, 1);
1059
1060     canvas->cd(1);
1061     InitPad();
1062     posX->DrawCopy();
1063
1064     canvas->cd(2);
1065     InitPad();
1066     posY->DrawCopy();
1067
1068     canvas->cd(3);
1069     InitPad();
1070     posZ->DrawCopy();
1071
1072     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1073     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1074   }
1075 }
1076
1077 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
1078 {
1079   TFile::Open(fileName);
1080   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1081   dNdEtaCorrection->LoadHistograms();
1082
1083   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1084   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1085
1086   gene->GetZaxis()->SetRangeUser(0.3, 10);
1087   meas->GetZaxis()->SetRangeUser(0.3, 10);
1088   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
1089   gene->GetZaxis()->SetRange(0, 0);
1090   meas->GetZaxis()->SetRange(0, 0);
1091
1092   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1093   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1094   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
1095   gene->GetYaxis()->SetRange(0, 0);
1096   meas->GetYaxis()->SetRange(0, 0);
1097
1098   gene->GetXaxis()->SetRangeUser(-10, 10);
1099   meas->GetXaxis()->SetRangeUser(-10, 10);
1100   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
1101   gene->GetXaxis()->SetRange(0, 0);
1102   meas->GetXaxis()->SetRange(0, 0);
1103 }
1104
1105 void Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1106 {
1107   gSystem->Load("libPWG0base");
1108
1109   Track2Particle2DCreatePlots(fileName);
1110
1111   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1112   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1113   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
1114
1115   /* this reads them from the file
1116   TH2* corrYX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx"));
1117   TH2* corrZX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx"));
1118   TH2* corrZY = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy"));*/
1119
1120   Prepare2DPlot(corrYX);
1121   Prepare2DPlot(corrZX);
1122   Prepare2DPlot(corrZY);
1123
1124   const char* title = "";
1125   corrYX->SetTitle(title);
1126   corrZX->SetTitle(title);
1127   corrZY->SetTitle(title);
1128
1129   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1130   canvas->Divide(3, 1);
1131
1132   canvas->cd(1);
1133   InitPadCOLZ();
1134   corrYX->Draw("COLZ");
1135
1136   canvas->cd(2);
1137   InitPadCOLZ();
1138   corrZX->Draw("COLZ");
1139
1140   canvas->cd(3);
1141   InitPadCOLZ();
1142   corrZY->Draw("COLZ");
1143
1144   canvas->SaveAs(Form("Track2Particle2D_%d.gif", gMax));
1145   canvas->SaveAs(Form("Track2Particle2D_%d.eps", gMax));
1146 }
1147
1148 void CompareTrack2Particle2D()
1149 {
1150   gSystem->Load("libPWG0base");
1151
1152   Track2Particle2DCreatePlots("correction_maponly-positive.root");
1153
1154   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("pos_yx"));
1155   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("pos_zx"));
1156   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("pos_zy"));
1157
1158   Track2Particle2DCreatePlots("correction_maponly-negative.root");
1159
1160   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("neg_yx"));
1161   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("neg_zx"));
1162   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("neg_zy"));
1163
1164   posYX->Divide(negYX);
1165   posZX->Divide(negZX);
1166   posZY->Divide(negZY);
1167
1168   Prepare2DPlot(posYX);
1169   Prepare2DPlot(posZX);
1170   Prepare2DPlot(posZY);
1171
1172   Float_t min = 0.8;
1173   Float_t max = 1.2;
1174
1175   posYX->SetMinimum(min);
1176   posYX->SetMaximum(max);
1177   posZX->SetMinimum(min);
1178   posZX->SetMaximum(max);
1179   posZY->SetMinimum(min);
1180   posZY->SetMaximum(max);
1181
1182   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1183   canvas->Divide(3, 1);
1184
1185   canvas->cd(1);
1186   InitPadCOLZ();
1187   posYX->Draw("COLZ");
1188
1189   canvas->cd(2);
1190   InitPadCOLZ();
1191   posZX->Draw("COLZ");
1192
1193   canvas->cd(3);
1194   InitPadCOLZ();
1195   posZY->Draw("COLZ");
1196
1197   canvas->SaveAs("CompareTrack2Particle2D.gif");
1198   canvas->SaveAs("CompareTrack2Particle2D.eps");
1199 }
1200
1201 void Track2Particle3D()
1202 {
1203   // get left margin proper
1204
1205   TFile* file = TFile::Open("correction_map.root");
1206
1207   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1208
1209   corr->SetTitle("Correction Factor");
1210   SetRanges(corr->GetZaxis());
1211
1212   Prepare3DPlot(corr);
1213
1214   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1215   canvas->SetTheta(29.428);
1216   canvas->SetPhi(16.5726);
1217
1218   corr->Draw();
1219
1220   canvas->SaveAs("Track2Particle3D.gif");
1221   canvas->SaveAs("Track2Particle3D.eps");
1222 }
1223
1224 void Track2Particle3DAll()
1225 {
1226   TFile* file = TFile::Open("correction_map.root");
1227
1228   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1229   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1230   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1231
1232   gene->SetTitle("Generated Particles");
1233   meas->SetTitle("Measured Tracks");
1234   corr->SetTitle("Correction Factor");
1235
1236   Prepare3DPlot(gene);
1237   Prepare3DPlot(meas);
1238   Prepare3DPlot(corr);
1239
1240   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1241   canvas->Divide(3, 1);
1242
1243   canvas->cd(1);
1244   InitPad();
1245   gene->Draw();
1246
1247   canvas->cd(2);
1248   meas->Draw();
1249
1250   canvas->cd(3);
1251   corr->Draw();
1252
1253   canvas->SaveAs("Track2Particle3DAll.gif");
1254   canvas->SaveAs("Track2Particle3DAll.eps");
1255 }
1256
1257 void MultiplicityMC(Int_t xRangeMax = 50)
1258 {
1259   TFile* file = TFile::Open("multiplicityMC.root");
1260
1261   if (!file)
1262   {
1263     printf("multiplicityMC.root could not be opened.\n");
1264     return;
1265   }
1266
1267   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1268   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1269   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1270
1271   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1272   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1273   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1274   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1275   {
1276     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1277     proj->Fit("gaus", "0");
1278     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1279     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1280
1281     continue;
1282
1283     // draw for debugging
1284     new TCanvas;
1285     proj->DrawCopy();
1286     proj->GetFunction("gaus")->DrawCopy("SAME");
1287   }
1288
1289   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1290
1291   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1292   {
1293     Float_t mean = correction->GetBinContent(i);
1294     Float_t width = correctionWidth->GetBinContent(i);
1295
1296     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1297     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1298     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1299
1300     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1301     {
1302       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1303     }
1304   }
1305
1306   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1307   fMultiplicityESDCorrectedRebinned->Rebin(10);
1308   fMultiplicityESDCorrectedRebinned->Scale(0.1);
1309
1310   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1311   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1312   ratio->Divide(fMultiplicityMC);
1313
1314   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1315   ratio2->Divide(fMultiplicityMC);
1316
1317   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1318   canvas->Divide(3, 2);
1319
1320   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1321   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1322   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1323   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1324   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1325   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1326   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1327
1328   canvas->cd(1); //InitPad();
1329   fMultiplicityESD->Draw();
1330   fMultiplicityMC->SetLineColor(2);
1331   fMultiplicityMC->Draw("SAME");
1332
1333   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1334   legend->AddEntry(fMultiplicityESD, "ESD");
1335   legend->AddEntry(fMultiplicityMC, "MC");
1336   legend->Draw();
1337
1338   canvas->cd(2);
1339   fCorrelation->Draw("COLZ");
1340
1341   canvas->cd(3);
1342   correction->Draw();
1343   //correction->Fit("pol1");
1344   correctionWidth->SetLineColor(2);
1345   correctionWidth->Draw("SAME");
1346
1347   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1348   legend->AddEntry(correction, "#bar{x}");
1349   legend->AddEntry(correctionWidth, "#sigma");
1350   legend->Draw();
1351
1352   canvas->cd(4);
1353   ratio->Draw();
1354
1355   ratio2->SetLineColor(2);
1356   ratio2->Draw("SAME");
1357
1358   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1359   legend->AddEntry(ratio, "uncorrected");
1360   legend->AddEntry(ratio2, "corrected");
1361   legend->Draw();
1362
1363   canvas->cd(5);
1364   fMultiplicityESDCorrected->SetLineColor(kBlue);
1365   fMultiplicityESDCorrected->Draw();
1366   fMultiplicityMC->Draw("SAME");
1367   fMultiplicityESD->Draw("SAME");
1368
1369   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1370   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1371   legend->AddEntry(fMultiplicityMC, "MC");
1372   legend->AddEntry(fMultiplicityESD, "ESD");
1373   legend->Draw();
1374
1375   canvas->cd(6);
1376   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
1377   fMultiplicityESDCorrectedRebinned->Draw();
1378   fMultiplicityMC->Draw("SAME");
1379
1380   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1381   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1382   legend->AddEntry(fMultiplicityMC, "MC");
1383   legend->Draw();
1384
1385   canvas->SaveAs("MultiplicityMC.gif");
1386 }
1387
1388 void MultiplicityESD()
1389 {
1390   TFile* file = TFile::Open("multiplicityESD.root");
1391
1392   if (!file)
1393   {
1394     printf("multiplicityESD.root could not be opened.\n");
1395     return;
1396   }
1397
1398   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1399
1400   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1401
1402   fMultiplicityESD->Draw();
1403 }
1404
1405 void drawPlots(Int_t max)
1406 {
1407   gMax = max;
1408
1409   ptCutoff();
1410   TriggerBias();
1411   VtxRecon();
1412   Track2Particle2D();
1413   Track2Particle3D();
1414   ptSpectrum();
1415   dNdEta();
1416 }
1417
1418 void drawPlots()
1419 {
1420   drawPlots(5);
1421   drawPlots(2);
1422 }