a3698971b5a066728649d66f6fb9342e524000d7
[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 dNdEta(Bool_t onlyESD = kFALSE)
98 {
99   TFile* file = TFile::Open("analysis_esd.root");
100   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
101   TH1* histESDNoPt = (TH1*) file->Get("dndeta/dndeta_dNdEta_2");
102   TH1* histESDMB = (TH1*) file->Get("dndeta_mb/dndeta_mb_dNdEta_corrected_2");
103   TH1* histESDMBVtx = (TH1*) file->Get("dndeta_mbvtx/dndeta_mbvtx_dNdEta_corrected_2");
104
105   TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
106
107   Prepare1DPlot(histESD);
108   Prepare1DPlot(histESDNoPt);
109   Prepare1DPlot(histESDMB);
110   Prepare1DPlot(histESDMBVtx);
111
112   histESD->SetLineColor(0);
113   histESDMB->SetLineColor(0);
114   histESDMBVtx->SetLineColor(0);
115
116   histESD->SetMarkerColor(kRed);
117   histESDMB->SetMarkerColor(kBlue);
118   histESDMBVtx->SetMarkerColor(103);
119
120   histESD->SetMarkerStyle(20);
121   histESDMB->SetMarkerStyle(21);
122   histESDMBVtx->SetMarkerStyle(22);
123
124   TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 100, 0, histESDMBVtx->GetMaximum() * 1.1);
125   Prepare1DPlot(dummy);
126   dummy->SetStats(kFALSE);
127   dummy->SetXTitle("#eta");
128   dummy->SetYTitle("dN_{ch}/d#eta");
129   dummy->GetYaxis()->SetTitleOffset(1);
130
131   histESDMBVtx->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
132   histESDMB->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
133   histESD->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
134   histESDNoPt->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
135
136   dummy->DrawCopy();
137   histESDMBVtx->Draw("SAME");
138   histESDMB->Draw("SAME");
139   histESD->Draw("SAME");
140
141   canvas->SaveAs("dNdEta1.gif");
142   canvas->SaveAs("dNdEta1.eps");
143
144   if (onlyESD)
145     return;
146
147   TFile* file2 = TFile::Open("analysis_mc.root");
148   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_dNdEta_corrected_2")->Clone("cloned");
149
150   gSystem->Load("libPWG0base");
151   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
152   fdNdEtaAnalysis->LoadHistograms();
153   fdNdEtaAnalysis->Finish(0, 0.3);
154   TH1* histMCPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(2);
155
156   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
157
158   Prepare1DPlot(histMC);
159   Prepare1DPlot(histMCPtCut);
160
161   histMC->SetLineColor(kBlue);
162   histMCPtCut->SetLineColor(104);
163   histESDNoPt->SetLineColor(102);
164
165   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
166   Prepare1DPlot(dummy2);
167   dummy2->GetYaxis()->SetRangeUser(0, histESD->GetMaximum() * 1.1);
168
169   dummy2->DrawCopy();
170   histMC->Draw("SAME");
171 //  histMC->Draw();
172   histESD->Draw("SAME");
173   histESDNoPt->Draw("SAME");
174   histMCPtCut->Draw("SAME");
175
176   canvas2->SaveAs("dNdEta2.gif");
177   canvas2->SaveAs("dNdEta2.eps");
178
179   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 500);
180   //InitPad();
181   gPad->SetBottomMargin(0.12);
182
183   TLegend* legend = new TLegend(0.35, 0.2, 0.6, 0.4);
184   legend->SetFillColor(0);
185   legend->AddEntry(histESDMBVtx, "triggered, vertex");
186   legend->AddEntry(histESDMB, "triggered");
187   legend->AddEntry(histESD, "all events");
188   legend->AddEntry(histMC, "MC prediction");
189
190   dummy->DrawCopy();
191   histESDMBVtx->Draw("SAME");
192   histESDMB->Draw("SAME");
193   histESD->Draw("SAME");
194   histMC->Draw("SAME");
195
196   legend->Draw();
197
198   canvas3->SaveAs("dNdEta.gif");
199   canvas3->SaveAs("dNdEta.eps");
200
201   TH1* ratio = (TH1*) histMC->Clone("ratio");
202   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
203
204   ratio->Divide(histESD);
205   ratioNoPt->Divide(histESDNoPt);
206
207   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
208
209   ratio->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
210
211   ratio->SetLineColor(1);
212   ratioNoPt->SetLineColor(2);
213
214   ratio->Draw();
215   ratioNoPt->Draw("SAME");
216
217 }
218
219 void ptSpectrum()
220 {
221   TFile* file = TFile::Open("analysis_esd.root");
222   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
223
224   TFile* file2 = TFile::Open("analysis_mc.root");
225   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
226
227   TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
228   InitPad();
229   gPad->SetLogy();
230
231   Prepare1DPlot(histMC);
232   Prepare1DPlot(histESD);
233
234   histESD->SetTitle("");
235   histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
236   histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
237
238   histMC->SetLineColor(kBlue);
239   histESD->SetLineColor(kRed);
240
241   histESD->GetYaxis()->SetTitleOffset(1.5);
242   histESD->GetXaxis()->SetRangeUser(0, 4.9999);
243
244   histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
245
246   histESD->Draw();
247   histMC->Draw("SAME");
248
249   canvas->SaveAs("ptSpectrum.gif");
250   canvas->SaveAs("ptSpectrum.eps");
251 }
252
253 void ptCutoff()
254 {
255   gSystem->Load("libPWG0base");
256
257   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
258   dNdEtaCorrection->LoadHistograms("correction_map.root","dndeta_correction");
259
260   dNdEtaCorrection->GetMeasuredFraction(0.3, -100, kTRUE);
261
262   TH1* hist = dynamic_cast<TH1*> (gROOT->FindObject("gene_dndeta_correction_nTrackToNPart_pt")->Clone("ptcutoff"));
263
264   hist->GetXaxis()->SetRangeUser(0, 0.9999);
265   hist->SetMinimum(0);
266
267   hist->SetTitle("Generated Particles");
268   Prepare1DPlot(hist);
269
270   TCanvas* canvas = new TCanvas("ptCutoff", "ptCutoff", 700, 500);
271   hist->DrawCopy();
272
273   TLine* line = new TLine(0.3, 0 - hist->GetMaximum() * 0, 0.3, hist->GetMaximum() * 1.1);
274   line->SetLineWidth(3);
275   line->SetLineColor(kRed);
276   line->Draw();
277
278   canvas->SaveAs("ptCutoff.gif");
279   canvas->SaveAs("ptCutoff.eps");
280
281   TH1F* factor = new TH1F("factor", ";#eta;correction factor", 10, -1, 1.000001);
282   for (Float_t eta = -0.9; eta<1; eta += 0.2)
283     factor->Fill(eta, 1.0 / dNdEtaCorrection->GetMeasuredFraction(0.3, eta, kFALSE));
284
285   TCanvas* canvas = new TCanvas("ptCutoff_factor", "ptCutoff_factor", 700, 500);
286   InitPad();
287
288   Prepare1DPlot(factor);
289   factor->GetYaxis()->SetRangeUser(1, 2);
290   factor->GetYaxis()->SetTitleOffset(1);
291   factor->Draw();
292
293   canvas->SaveAs("ptCutoff_factor.eps");
294 }
295
296 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
297 {
298   TFile* file = TFile::Open(fileName);
299
300   TH2* corrTrigger = dynamic_cast<TH2*> (file->Get(Form("%s/corr_%s_trigger", folder, folder)));
301   TH2* corrVtx = dynamic_cast<TH2*> (file->Get(Form("%s/corr_%s_vtxReco", folder, folder)));
302
303   Prepare2DPlot(corrTrigger);
304   corrTrigger->SetTitle("a) Trigger bias correction");
305
306   Prepare2DPlot(corrVtx);
307   corrVtx->SetTitle("b) Vertex reconstruction correction");
308
309   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
310   canvas->Divide(2, 1);
311
312   canvas->cd(1);
313   InitPadCOLZ();
314   corrTrigger->DrawCopy("COLZ");
315
316   canvas->cd(2);
317   InitPadCOLZ();
318   corrVtx->DrawCopy("COLZ");
319
320   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
321   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
322
323   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
324   canvas->Divide(2, 1);
325
326   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
327   corrVtx->GetYaxis()->SetRangeUser(0, 5);
328
329   canvas->cd(1);
330   InitPadCOLZ();
331   corrTrigger->DrawCopy("COLZ");
332
333   canvas->cd(2);
334   InitPadCOLZ();
335   corrVtx->DrawCopy("COLZ");
336
337   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
338   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
339 }
340
341 void TriggerBias(const char* fileName = "correction_map.root")
342 {
343   TFile* file = TFile::Open(fileName);
344
345   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_trigger"));
346
347   Prepare2DPlot(corr);
348   corr->SetTitle("Trigger bias correction");
349
350   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
351   InitPadCOLZ();
352   corr->DrawCopy("COLZ");
353
354   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
355   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
356
357   corr->GetYaxis()->SetRangeUser(0, 5);
358
359   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
360   InitPadCOLZ();
361   corr->DrawCopy("COLZ");
362
363   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
364   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
365 }
366
367 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
368 {
369   gSystem->Load("libPWG0base");
370
371   TFile* file = TFile::Open(fileName);
372   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
373   dNdEtaCorrection->LoadHistograms(fileName, folderName);
374
375   TH1* hist = dNdEtaCorrection->GetTriggerCorrection()->Get1DCorrection("x");
376
377   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 500, 500);
378   InitPad();
379
380   Prepare1DPlot(hist);
381   hist->SetTitle("");
382   hist->GetYaxis()->SetTitle("correction factor");
383   hist->GetYaxis()->SetRangeUser(1, 1.5);
384   hist->GetYaxis()->SetTitleOffset(1.6);
385   hist->Draw();
386
387   canvas->SaveAs("TriggerBias1D.eps");
388 }
389
390 void VtxRecon()
391 {
392   TFile* file = TFile::Open("correction_map.root");
393
394   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
395
396   Prepare2DPlot(corr);
397   corr->SetTitle("Vertex reconstruction correction");
398
399   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
400   InitPadCOLZ();
401   corr->DrawCopy("COLZ");
402
403   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
404   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
405
406   corr->GetYaxis()->SetRangeUser(0, 5);
407
408   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
409   InitPadCOLZ();
410   corr->DrawCopy("COLZ");
411
412   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
413   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
414 }
415
416 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
417 {
418   gSystem->Load("libPWG0base");
419
420   TFile* file = TFile::Open(fileName);
421   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
422   dNdEtaCorrection->LoadHistograms(fileName, folderName);
423
424   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->Get1DCorrection("x");
425
426   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 500, 500);
427   InitPad();
428
429   Prepare1DPlot(hist);
430   hist->SetTitle("");
431   hist->GetYaxis()->SetTitle("correction factor");
432   hist->GetYaxis()->SetRangeUser(1, 1.8);
433   hist->GetYaxis()->SetTitleOffset(1.6);
434   hist->Draw();
435
436   canvas->SaveAs("VtxRecon1D.eps");
437 }
438
439 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
440 {
441   gSystem->Load("libPWG0base");
442
443   TFile::Open(fileName);
444   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
445   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
446
447   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
448   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
449
450   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
451   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
452
453   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
454   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
455   gene->GetXaxis()->SetRangeUser(-10, 10);
456   meas->GetXaxis()->SetRangeUser(-10, 10);
457
458   Float_t eff1 = gene->Integral() / meas->Integral();
459   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
460
461   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
462
463   gene->GetZaxis()->SetRangeUser(0.3, 10);
464   meas->GetZaxis()->SetRangeUser(0.3, 10);
465
466   Float_t eff2 = gene->Integral() / meas->Integral();
467   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
468
469   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
470
471   gene->GetZaxis()->SetRangeUser(0.3, 1);
472   meas->GetZaxis()->SetRangeUser(0.3, 1);
473
474   Float_t eff3 = gene->Integral() / meas->Integral();
475   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
476
477   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
478 }
479
480 void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 10)
481 {
482   TFile::Open(fileName);
483   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
484   dNdEtaCorrection->LoadHistograms(fileName, folderName);
485
486   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
487   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
488
489   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
490   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
491   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
492   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
493   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kTRUE);
494   gene->GetYaxis()->SetRange(0, 0);
495   meas->GetYaxis()->SetRange(0, 0);
496
497   gene->GetXaxis()->SetRangeUser(-10, 10);
498   meas->GetXaxis()->SetRangeUser(-10, 10);
499   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kTRUE);
500   gene->GetZaxis()->SetRange(0, 0);
501   meas->GetZaxis()->SetRange(0, 0);
502
503   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
504   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
505   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kTRUE);
506 }
507
508 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
509 {
510   gSystem->Load("libPWG0base");
511
512   Track2Particle1DCreatePlots(fileName, folder, upperPtLimit);
513
514   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", folder, folder)));
515   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", folder, folder)));
516   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", folder, folder)));
517
518   Prepare1DPlot(corrX);
519   Prepare1DPlot(corrY);
520   Prepare1DPlot(corrZ);
521
522   //corrX->SetTitle("a) z projection");
523   corrY->SetTitle("a) #eta projection");
524   corrZ->SetTitle("b) p_{T} projection");
525
526   corrY->GetYaxis()->SetTitle("correction factor");
527   corrZ->GetYaxis()->SetTitle("correction factor");
528
529   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
530
531   TString canvasName;
532   canvasName.Form("Track2Particle1D_%s", folder);
533   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
534   canvas->Divide(3, 1);
535
536   canvas->cd(1);
537   InitPad();
538   corrX->DrawCopy();
539
540   canvas->cd(2);
541   InitPad();
542   corrY->DrawCopy();
543
544   canvas->cd(3);
545   InitPad();
546   corrZ->DrawCopy();
547
548   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
549   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
550
551   canvasName.Form("Track2Particle1D_%s_etapt", folder);
552   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
553   canvas->Divide(2, 1);
554
555   canvas->cd(1);
556   InitPad();
557   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
558   corrY->GetYaxis()->SetRangeUser(1, 1.5);
559   corrY->GetYaxis()->SetTitleOffset(1.5);
560   corrY->DrawCopy();
561
562   canvas->cd(2);
563   InitPad();
564   corrZ->GetYaxis()->SetRangeUser(1, 1.5);
565   corrZ->GetYaxis()->SetTitleOffset(1.5);
566   corrZ->DrawCopy();
567
568   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
569 }
570
571 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
572 {
573   gSystem->Load("libPWG0base");
574
575   Track2Particle1DCreatePlots("correction_maponly-positive.root", "dndeta_correction", upperPtLimit);
576
577   TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_x_div_meas_nTrackToNPart_x")->Clone("pos_x"));
578   TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_y_div_meas_nTrackToNPart_y")->Clone("pos_y"));
579   TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_z_div_meas_nTrackToNPart_z")->Clone("pos_z"));
580
581   Track2Particle1DCreatePlots("correction_maponly-negative.root", "dndeta_correction", upperPtLimit);
582
583   TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_x_div_meas_nTrackToNPart_x")->Clone("neg_x"));
584   TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_y_div_meas_nTrackToNPart_y")->Clone("neg_y"));
585   TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_z_div_meas_nTrackToNPart_z")->Clone("neg_z"));
586
587   //printf("%f %f %f %f\n", posX->GetBinContent(20), posX->GetBinError(20), negX->GetBinContent(20), negX->GetBinError(20));
588
589   posX->Divide(negX);
590   posY->Divide(negY);
591   posZ->Divide(negZ);
592
593   //printf("%f %f\n", posX->GetBinContent(20), posX->GetBinError(20));
594
595   Prepare1DPlot(posX);
596   Prepare1DPlot(posY);
597   Prepare1DPlot(posZ);
598
599   Float_t min = 0.8;
600   Float_t max = 1.2;
601
602   posX->SetMinimum(min);
603   posX->SetMaximum(max);
604   posY->SetMinimum(min);
605   posY->SetMaximum(max);
606   posZ->SetMinimum(min);
607   posZ->SetMaximum(max);
608
609   posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
610
611   posX->GetYaxis()->SetTitleOffset(1.7);
612   posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
613   posY->GetYaxis()->SetTitleOffset(1.7);
614   posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
615   posZ->GetYaxis()->SetTitleOffset(1.7);
616   posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
617
618   TCanvas* canvas = new TCanvas("CompareTrack2Particle1D", "CompareTrack2Particle1D", 1200, 400);
619   canvas->Divide(3, 1);
620
621   canvas->cd(1);
622   InitPad();
623   posX->Draw();
624
625   canvas->cd(2);
626   InitPad();
627   posY->Draw();
628
629   canvas->cd(3);
630   InitPad();
631   posZ->Draw();
632
633   canvas->SaveAs(Form("CompareTrack2Particle1D_%f.gif", upperPtLimit));
634   canvas->SaveAs(Form("CompareTrack2Particle1D_%f.eps", upperPtLimit));
635 }
636
637 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
638 {
639   TFile::Open(fileName);
640   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
641   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
642
643   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
644   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
645
646   gene->GetZaxis()->SetRangeUser(0.3, 10);
647   meas->GetZaxis()->SetRangeUser(0.3, 10);
648   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
649   gene->GetZaxis()->SetRange(0, 0);
650   meas->GetZaxis()->SetRange(0, 0);
651
652   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
653   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
654   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
655   gene->GetYaxis()->SetRange(0, 0);
656   meas->GetYaxis()->SetRange(0, 0);
657
658   gene->GetXaxis()->SetRangeUser(-10, 10);
659   meas->GetXaxis()->SetRangeUser(-10, 10);
660   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
661   gene->GetXaxis()->SetRange(0, 0);
662   meas->GetXaxis()->SetRange(0, 0);
663 }
664
665 void Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
666 {
667   gSystem->Load("libPWG0base");
668
669   Track2Particle2DCreatePlots(fileName);
670
671   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_yx_div_meas_%s_nTrackToNPart_yx", folder, folder)));
672   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_zx_div_meas_%s_nTrackToNPart_zx", folder, folder)));
673   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_zy_div_meas_%s_nTrackToNPart_zy", folder, folder)));
674
675   /* this reads them from the file
676   TH2* corrYX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx"));
677   TH2* corrZX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx"));
678   TH2* corrZY = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy"));*/
679
680   Prepare2DPlot(corrYX);
681   Prepare2DPlot(corrZX);
682   Prepare2DPlot(corrZY);
683
684   const char* title = "";
685   corrYX->SetTitle(title);
686   corrZX->SetTitle(title);
687   corrZY->SetTitle(title);
688
689   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
690   canvas->Divide(3, 1);
691
692   canvas->cd(1);
693   InitPadCOLZ();
694   corrYX->Draw("COLZ");
695
696   canvas->cd(2);
697   InitPadCOLZ();
698   corrZX->Draw("COLZ");
699
700   canvas->cd(3);
701   InitPadCOLZ();
702   corrZY->Draw("COLZ");
703
704   canvas->SaveAs(Form("Track2Particle2D_%d.gif", gMax));
705   canvas->SaveAs(Form("Track2Particle2D_%d.eps", gMax));
706 }
707
708 void CompareTrack2Particle2D()
709 {
710   gSystem->Load("libPWG0base");
711
712   Track2Particle2DCreatePlots("correction_maponly-positive.root");
713
714   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("pos_yx"));
715   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("pos_zx"));
716   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("pos_zy"));
717
718   Track2Particle2DCreatePlots("correction_maponly-negative.root");
719
720   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("neg_yx"));
721   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("neg_zx"));
722   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("neg_zy"));
723
724   posYX->Divide(negYX);
725   posZX->Divide(negZX);
726   posZY->Divide(negZY);
727
728   Prepare2DPlot(posYX);
729   Prepare2DPlot(posZX);
730   Prepare2DPlot(posZY);
731
732   Float_t min = 0.8;
733   Float_t max = 1.2;
734
735   posYX->SetMinimum(min);
736   posYX->SetMaximum(max);
737   posZX->SetMinimum(min);
738   posZX->SetMaximum(max);
739   posZY->SetMinimum(min);
740   posZY->SetMaximum(max);
741
742   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
743   canvas->Divide(3, 1);
744
745   canvas->cd(1);
746   InitPadCOLZ();
747   posYX->Draw("COLZ");
748
749   canvas->cd(2);
750   InitPadCOLZ();
751   posZX->Draw("COLZ");
752
753   canvas->cd(3);
754   InitPadCOLZ();
755   posZY->Draw("COLZ");
756
757   canvas->SaveAs("CompareTrack2Particle2D.gif");
758   canvas->SaveAs("CompareTrack2Particle2D.eps");
759 }
760
761 void Track2Particle3D()
762 {
763   // get left margin proper
764
765   TFile* file = TFile::Open("correction_map.root");
766
767   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
768
769   corr->SetTitle("Correction Factor");
770   SetRanges(corr->GetZaxis());
771
772   Prepare3DPlot(corr);
773
774   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
775   canvas->SetTheta(29.428);
776   canvas->SetPhi(16.5726);
777
778   corr->Draw();
779
780   canvas->SaveAs("Track2Particle3D.gif");
781   canvas->SaveAs("Track2Particle3D.eps");
782 }
783
784 void Track2Particle3DAll()
785 {
786   TFile* file = TFile::Open("correction_map.root");
787
788   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
789   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
790   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
791
792   gene->SetTitle("Generated Particles");
793   meas->SetTitle("Measured Tracks");
794   corr->SetTitle("Correction Factor");
795
796   Prepare3DPlot(gene);
797   Prepare3DPlot(meas);
798   Prepare3DPlot(corr);
799
800   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
801   canvas->Divide(3, 1);
802
803   canvas->cd(1);
804   InitPad();
805   gene->Draw();
806
807   canvas->cd(2);
808   meas->Draw();
809
810   canvas->cd(3);
811   corr->Draw();
812
813   canvas->SaveAs("Track2Particle3DAll.gif");
814   canvas->SaveAs("Track2Particle3DAll.eps");
815 }
816
817 void MultiplicityMC(Int_t xRangeMax = 50)
818 {
819   TFile* file = TFile::Open("multiplicityMC.root");
820
821   if (!file)
822   {
823     printf("multiplicityMC.root could not be opened.\n");
824     return;
825   }
826
827   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
828   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
829   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
830
831   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
832   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
833   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
834   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
835   {
836     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
837     proj->Fit("gaus", "0");
838     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
839     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
840
841     continue;
842
843     // draw for debugging
844     new TCanvas;
845     proj->DrawCopy();
846     proj->GetFunction("gaus")->DrawCopy("SAME");
847   }
848
849   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
850
851   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
852   {
853     Float_t mean = correction->GetBinContent(i);
854     Float_t width = correctionWidth->GetBinContent(i);
855
856     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
857     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
858     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
859
860     for (Int_t j=fillBegin; j <= fillEnd; ++j)
861     {
862       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
863     }
864   }
865
866   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
867   fMultiplicityESDCorrectedRebinned->Rebin(10);
868   fMultiplicityESDCorrectedRebinned->Scale(0.1);
869
870   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
871   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
872   ratio->Divide(fMultiplicityMC);
873
874   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
875   ratio2->Divide(fMultiplicityMC);
876
877   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
878   canvas->Divide(3, 2);
879
880   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
881   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
882   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
883   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
884   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
885   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
886   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
887
888   canvas->cd(1); //InitPad();
889   fMultiplicityESD->Draw();
890   fMultiplicityMC->SetLineColor(2);
891   fMultiplicityMC->Draw("SAME");
892
893   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
894   legend->AddEntry(fMultiplicityESD, "ESD");
895   legend->AddEntry(fMultiplicityMC, "MC");
896   legend->Draw();
897
898   canvas->cd(2);
899   fCorrelation->Draw("COLZ");
900
901   canvas->cd(3);
902   correction->Draw();
903   //correction->Fit("pol1");
904   correctionWidth->SetLineColor(2);
905   correctionWidth->Draw("SAME");
906
907   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
908   legend->AddEntry(correction, "#bar{x}");
909   legend->AddEntry(correctionWidth, "#sigma");
910   legend->Draw();
911
912   canvas->cd(4);
913   ratio->Draw();
914
915   ratio2->SetLineColor(2);
916   ratio2->Draw("SAME");
917
918   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
919   legend->AddEntry(ratio, "uncorrected");
920   legend->AddEntry(ratio2, "corrected");
921   legend->Draw();
922
923   canvas->cd(5);
924   fMultiplicityESDCorrected->SetLineColor(kBlue);
925   fMultiplicityESDCorrected->Draw();
926   fMultiplicityMC->Draw("SAME");
927   fMultiplicityESD->Draw("SAME");
928
929   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
930   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
931   legend->AddEntry(fMultiplicityMC, "MC");
932   legend->AddEntry(fMultiplicityESD, "ESD");
933   legend->Draw();
934
935   canvas->cd(6);
936   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
937   fMultiplicityESDCorrectedRebinned->Draw();
938   fMultiplicityMC->Draw("SAME");
939
940   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
941   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
942   legend->AddEntry(fMultiplicityMC, "MC");
943   legend->Draw();
944
945   canvas->SaveAs("MultiplicityMC.gif");
946 }
947
948 void MultiplicityESD()
949 {
950   TFile* file = TFile::Open("multiplicityESD.root");
951
952   if (!file)
953   {
954     printf("multiplicityESD.root could not be opened.\n");
955     return;
956   }
957
958   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
959
960   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
961
962   fMultiplicityESD->Draw();
963 }
964
965 void drawPlots(Int_t max)
966 {
967   gMax = max;
968
969   ptCutoff();
970   TriggerBias();
971   VtxRecon();
972   Track2Particle2D();
973   Track2Particle3D();
974   ptSpectrum();
975   dNdEta();
976 }
977
978 void drawPlots()
979 {
980   drawPlots(5);
981   drawPlots(2);
982 }