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