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