updating makefile to include evgen headers to compile systematics selector
[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   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
310   corrVtx->GetYaxis()->SetTitle("Multiplicity");
311
312   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
313   canvas->Divide(2, 1);
314
315   canvas->cd(1);
316   InitPadCOLZ();
317   corrTrigger->DrawCopy("COLZ");
318
319   canvas->cd(2);
320   InitPadCOLZ();
321   corrVtx->DrawCopy("COLZ");
322
323   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
324   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
325
326   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
327   canvas->Divide(2, 1);
328
329   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
330   corrVtx->GetYaxis()->SetRangeUser(0, 5);
331
332   canvas->cd(1);
333   InitPadCOLZ();
334   corrTrigger->DrawCopy("COLZ");
335
336   canvas->cd(2);
337   InitPadCOLZ();
338   corrVtx->DrawCopy("COLZ");
339
340   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
341   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
342 }
343
344 void TriggerBias(const char* fileName = "correction_map.root")
345 {
346   TFile* file = TFile::Open(fileName);
347
348   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
349
350   Prepare2DPlot(corr);
351   corr->SetTitle("Trigger bias correction");
352
353   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
354   InitPadCOLZ();
355   corr->DrawCopy("COLZ");
356
357   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
358   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
359
360   corr->GetYaxis()->SetRangeUser(0, 5);
361
362   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
363   InitPadCOLZ();
364   corr->DrawCopy("COLZ");
365
366   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
367   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
368 }
369
370 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
371 {
372   gSystem->Load("libPWG0base");
373
374   TFile* file = TFile::Open(fileName);
375   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
376   dNdEtaCorrection->LoadHistograms(fileName, folderName);
377
378   TH1* hist = dNdEtaCorrection->GetTriggerCorrection()->Get1DCorrection("x");
379
380   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 500, 500);
381   InitPad();
382
383   Prepare1DPlot(hist);
384   hist->SetTitle("");
385   hist->GetYaxis()->SetTitle("correction factor");
386   hist->GetYaxis()->SetRangeUser(1, 1.5);
387   hist->GetYaxis()->SetTitleOffset(1.6);
388   hist->Draw();
389
390   canvas->SaveAs("TriggerBias1D.eps");
391 }
392
393 void VtxRecon()
394 {
395   TFile* file = TFile::Open("correction_map.root");
396
397   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
398
399   Prepare2DPlot(corr);
400   corr->SetTitle("Vertex reconstruction correction");
401
402   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
403   InitPadCOLZ();
404   corr->DrawCopy("COLZ");
405
406   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
407   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
408
409   corr->GetYaxis()->SetRangeUser(0, 5);
410
411   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
412   InitPadCOLZ();
413   corr->DrawCopy("COLZ");
414
415   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
416   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
417 }
418
419 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
420 {
421   gSystem->Load("libPWG0base");
422
423   TFile* file = TFile::Open(fileName);
424   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
425   dNdEtaCorrection->LoadHistograms(fileName, folderName);
426
427   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->Get1DCorrection("x");
428
429   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 500, 500);
430   InitPad();
431
432   Prepare1DPlot(hist);
433   hist->SetTitle("");
434   hist->GetYaxis()->SetTitle("correction factor");
435   hist->GetYaxis()->SetRangeUser(1, 1.8);
436   hist->GetYaxis()->SetTitleOffset(1.6);
437   hist->Draw();
438
439   canvas->SaveAs("VtxRecon1D.eps");
440 }
441
442 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
443 {
444   gSystem->Load("libPWG0base");
445
446   TFile::Open(fileName);
447   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
448   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
449
450   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
451   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
452
453   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
454   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
455
456   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
457   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
458   gene->GetXaxis()->SetRangeUser(-10, 10);
459   meas->GetXaxis()->SetRangeUser(-10, 10);
460
461   Float_t eff1 = gene->Integral() / meas->Integral();
462   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
463
464   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
465
466   gene->GetZaxis()->SetRangeUser(0.3, 10);
467   meas->GetZaxis()->SetRangeUser(0.3, 10);
468
469   Float_t eff2 = gene->Integral() / meas->Integral();
470   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
471
472   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
473
474   gene->GetZaxis()->SetRangeUser(0.3, 1);
475   meas->GetZaxis()->SetRangeUser(0.3, 1);
476
477   Float_t eff3 = gene->Integral() / meas->Integral();
478   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
479
480   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
481 }
482
483 void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 10)
484 {
485   TFile::Open(fileName);
486   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
487   dNdEtaCorrection->LoadHistograms(fileName, folderName);
488
489   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
490   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
491
492   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
493   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
494   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
495   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
496   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kTRUE);
497   gene->GetYaxis()->SetRange(0, 0);
498   meas->GetYaxis()->SetRange(0, 0);
499
500   gene->GetXaxis()->SetRangeUser(-10, 10);
501   meas->GetXaxis()->SetRangeUser(-10, 10);
502   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kTRUE);
503   gene->GetZaxis()->SetRange(0, 0);
504   meas->GetZaxis()->SetRange(0, 0);
505
506   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
507   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
508   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kTRUE);
509 }
510
511 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
512 {
513   gSystem->Load("libPWG0base");
514
515   Track2Particle1DCreatePlots(fileName, folder, upperPtLimit);
516
517   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", folder, folder)));
518   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", folder, folder)));
519   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", folder, folder)));
520
521   Prepare1DPlot(corrX);
522   Prepare1DPlot(corrY);
523   Prepare1DPlot(corrZ);
524
525   //corrX->SetTitle("a) z projection");
526   corrY->SetTitle("a) #eta projection");
527   corrZ->SetTitle("b) p_{T} projection");
528
529   corrY->GetYaxis()->SetTitle("correction factor");
530   corrZ->GetYaxis()->SetTitle("correction factor");
531
532   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
533
534   TString canvasName;
535   canvasName.Form("Track2Particle1D_%s", folder);
536   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
537   canvas->Divide(3, 1);
538
539   canvas->cd(1);
540   InitPad();
541   corrX->DrawCopy();
542
543   canvas->cd(2);
544   InitPad();
545   corrY->DrawCopy();
546
547   canvas->cd(3);
548   InitPad();
549   corrZ->DrawCopy();
550
551   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
552   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
553
554   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
555
556   canvasName.Form("Track2Particle1D_%s_etapt", folder);
557   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
558   canvas->Divide(2, 1);
559
560   canvas->cd(1);
561   InitPad();
562   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
563   corrY->GetYaxis()->SetRangeUser(1, 1.5);
564   corrY->GetYaxis()->SetTitleOffset(1.5);
565   corrY->DrawCopy();
566   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
567   pave->AddText("|z| < 10 cm");
568   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
569   pave->Draw();
570
571   canvas->cd(2);
572   InitPad();
573   gPad->SetLogx();
574   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
575   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
576   corrZ->GetYaxis()->SetTitleOffset(1.5);
577   corrZ->DrawCopy();
578   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
579   pave->AddText("|z| < 10 cm");
580   pave->AddText("|#eta| < 0.8");
581   pave->Draw();
582
583   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
584 }
585
586 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
587 {
588   gSystem->Load("libPWG0base");
589
590   Track2Particle1DCreatePlots("correction_maponly-positive.root", "dndeta_correction", upperPtLimit);
591
592   TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_x_div_meas_nTrackToNPart_x")->Clone("pos_x"));
593   TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_y_div_meas_nTrackToNPart_y")->Clone("pos_y"));
594   TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_z_div_meas_nTrackToNPart_z")->Clone("pos_z"));
595
596   Track2Particle1DCreatePlots("correction_maponly-negative.root", "dndeta_correction", upperPtLimit);
597
598   TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_x_div_meas_nTrackToNPart_x")->Clone("neg_x"));
599   TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_y_div_meas_nTrackToNPart_y")->Clone("neg_y"));
600   TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_z_div_meas_nTrackToNPart_z")->Clone("neg_z"));
601
602   //printf("%f %f %f %f\n", posX->GetBinContent(20), posX->GetBinError(20), negX->GetBinContent(20), negX->GetBinError(20));
603
604   posX->Divide(negX);
605   posY->Divide(negY);
606   posZ->Divide(negZ);
607
608   //printf("%f %f\n", posX->GetBinContent(20), posX->GetBinError(20));
609
610   Prepare1DPlot(posX);
611   Prepare1DPlot(posY);
612   Prepare1DPlot(posZ);
613
614   Float_t min = 0.8;
615   Float_t max = 1.2;
616
617   posX->SetMinimum(min);
618   posX->SetMaximum(max);
619   posY->SetMinimum(min);
620   posY->SetMaximum(max);
621   posZ->SetMinimum(min);
622   posZ->SetMaximum(max);
623
624   posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
625
626   posX->GetYaxis()->SetTitleOffset(1.7);
627   posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
628   posY->GetYaxis()->SetTitleOffset(1.7);
629   posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
630   posZ->GetYaxis()->SetTitleOffset(1.7);
631   posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
632
633   TCanvas* canvas = new TCanvas("CompareTrack2Particle1D", "CompareTrack2Particle1D", 1200, 400);
634   canvas->Divide(3, 1);
635
636   canvas->cd(1);
637   InitPad();
638   posX->Draw();
639
640   canvas->cd(2);
641   InitPad();
642   posY->Draw();
643
644   canvas->cd(3);
645   InitPad();
646   posZ->Draw();
647
648   canvas->SaveAs(Form("CompareTrack2Particle1D_%f.gif", upperPtLimit));
649   canvas->SaveAs(Form("CompareTrack2Particle1D_%f.eps", upperPtLimit));
650 }
651
652 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
653 {
654   TFile::Open(fileName);
655   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
656   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
657
658   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
659   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
660
661   gene->GetZaxis()->SetRangeUser(0.3, 10);
662   meas->GetZaxis()->SetRangeUser(0.3, 10);
663   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
664   gene->GetZaxis()->SetRange(0, 0);
665   meas->GetZaxis()->SetRange(0, 0);
666
667   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
668   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
669   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
670   gene->GetYaxis()->SetRange(0, 0);
671   meas->GetYaxis()->SetRange(0, 0);
672
673   gene->GetXaxis()->SetRangeUser(-10, 10);
674   meas->GetXaxis()->SetRangeUser(-10, 10);
675   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
676   gene->GetXaxis()->SetRange(0, 0);
677   meas->GetXaxis()->SetRange(0, 0);
678 }
679
680 void Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
681 {
682   gSystem->Load("libPWG0base");
683
684   Track2Particle2DCreatePlots(fileName);
685
686   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_yx_div_meas_%s_nTrackToNPart_yx", folder, folder)));
687   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_zx_div_meas_%s_nTrackToNPart_zx", folder, folder)));
688   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_zy_div_meas_%s_nTrackToNPart_zy", folder, folder)));
689
690   /* this reads them from the file
691   TH2* corrYX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx"));
692   TH2* corrZX = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx"));
693   TH2* corrZY = dynamic_cast<TH2*> (file->Get("dndeta_correction/gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy"));*/
694
695   Prepare2DPlot(corrYX);
696   Prepare2DPlot(corrZX);
697   Prepare2DPlot(corrZY);
698
699   const char* title = "";
700   corrYX->SetTitle(title);
701   corrZX->SetTitle(title);
702   corrZY->SetTitle(title);
703
704   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
705   canvas->Divide(3, 1);
706
707   canvas->cd(1);
708   InitPadCOLZ();
709   corrYX->Draw("COLZ");
710
711   canvas->cd(2);
712   InitPadCOLZ();
713   corrZX->Draw("COLZ");
714
715   canvas->cd(3);
716   InitPadCOLZ();
717   corrZY->Draw("COLZ");
718
719   canvas->SaveAs(Form("Track2Particle2D_%d.gif", gMax));
720   canvas->SaveAs(Form("Track2Particle2D_%d.eps", gMax));
721 }
722
723 void CompareTrack2Particle2D()
724 {
725   gSystem->Load("libPWG0base");
726
727   Track2Particle2DCreatePlots("correction_maponly-positive.root");
728
729   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("pos_yx"));
730   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("pos_zx"));
731   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("pos_zy"));
732
733   Track2Particle2DCreatePlots("correction_maponly-negative.root");
734
735   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("neg_yx"));
736   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("neg_zx"));
737   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("neg_zy"));
738
739   posYX->Divide(negYX);
740   posZX->Divide(negZX);
741   posZY->Divide(negZY);
742
743   Prepare2DPlot(posYX);
744   Prepare2DPlot(posZX);
745   Prepare2DPlot(posZY);
746
747   Float_t min = 0.8;
748   Float_t max = 1.2;
749
750   posYX->SetMinimum(min);
751   posYX->SetMaximum(max);
752   posZX->SetMinimum(min);
753   posZX->SetMaximum(max);
754   posZY->SetMinimum(min);
755   posZY->SetMaximum(max);
756
757   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
758   canvas->Divide(3, 1);
759
760   canvas->cd(1);
761   InitPadCOLZ();
762   posYX->Draw("COLZ");
763
764   canvas->cd(2);
765   InitPadCOLZ();
766   posZX->Draw("COLZ");
767
768   canvas->cd(3);
769   InitPadCOLZ();
770   posZY->Draw("COLZ");
771
772   canvas->SaveAs("CompareTrack2Particle2D.gif");
773   canvas->SaveAs("CompareTrack2Particle2D.eps");
774 }
775
776 void Track2Particle3D()
777 {
778   // get left margin proper
779
780   TFile* file = TFile::Open("correction_map.root");
781
782   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
783
784   corr->SetTitle("Correction Factor");
785   SetRanges(corr->GetZaxis());
786
787   Prepare3DPlot(corr);
788
789   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
790   canvas->SetTheta(29.428);
791   canvas->SetPhi(16.5726);
792
793   corr->Draw();
794
795   canvas->SaveAs("Track2Particle3D.gif");
796   canvas->SaveAs("Track2Particle3D.eps");
797 }
798
799 void Track2Particle3DAll()
800 {
801   TFile* file = TFile::Open("correction_map.root");
802
803   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
804   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
805   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
806
807   gene->SetTitle("Generated Particles");
808   meas->SetTitle("Measured Tracks");
809   corr->SetTitle("Correction Factor");
810
811   Prepare3DPlot(gene);
812   Prepare3DPlot(meas);
813   Prepare3DPlot(corr);
814
815   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
816   canvas->Divide(3, 1);
817
818   canvas->cd(1);
819   InitPad();
820   gene->Draw();
821
822   canvas->cd(2);
823   meas->Draw();
824
825   canvas->cd(3);
826   corr->Draw();
827
828   canvas->SaveAs("Track2Particle3DAll.gif");
829   canvas->SaveAs("Track2Particle3DAll.eps");
830 }
831
832 void MultiplicityMC(Int_t xRangeMax = 50)
833 {
834   TFile* file = TFile::Open("multiplicityMC.root");
835
836   if (!file)
837   {
838     printf("multiplicityMC.root could not be opened.\n");
839     return;
840   }
841
842   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
843   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
844   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
845
846   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
847   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
848   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
849   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
850   {
851     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
852     proj->Fit("gaus", "0");
853     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
854     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
855
856     continue;
857
858     // draw for debugging
859     new TCanvas;
860     proj->DrawCopy();
861     proj->GetFunction("gaus")->DrawCopy("SAME");
862   }
863
864   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
865
866   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
867   {
868     Float_t mean = correction->GetBinContent(i);
869     Float_t width = correctionWidth->GetBinContent(i);
870
871     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
872     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
873     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
874
875     for (Int_t j=fillBegin; j <= fillEnd; ++j)
876     {
877       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
878     }
879   }
880
881   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
882   fMultiplicityESDCorrectedRebinned->Rebin(10);
883   fMultiplicityESDCorrectedRebinned->Scale(0.1);
884
885   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
886   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
887   ratio->Divide(fMultiplicityMC);
888
889   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
890   ratio2->Divide(fMultiplicityMC);
891
892   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
893   canvas->Divide(3, 2);
894
895   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
896   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
897   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
898   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
899   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
900   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
901   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
902
903   canvas->cd(1); //InitPad();
904   fMultiplicityESD->Draw();
905   fMultiplicityMC->SetLineColor(2);
906   fMultiplicityMC->Draw("SAME");
907
908   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
909   legend->AddEntry(fMultiplicityESD, "ESD");
910   legend->AddEntry(fMultiplicityMC, "MC");
911   legend->Draw();
912
913   canvas->cd(2);
914   fCorrelation->Draw("COLZ");
915
916   canvas->cd(3);
917   correction->Draw();
918   //correction->Fit("pol1");
919   correctionWidth->SetLineColor(2);
920   correctionWidth->Draw("SAME");
921
922   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
923   legend->AddEntry(correction, "#bar{x}");
924   legend->AddEntry(correctionWidth, "#sigma");
925   legend->Draw();
926
927   canvas->cd(4);
928   ratio->Draw();
929
930   ratio2->SetLineColor(2);
931   ratio2->Draw("SAME");
932
933   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
934   legend->AddEntry(ratio, "uncorrected");
935   legend->AddEntry(ratio2, "corrected");
936   legend->Draw();
937
938   canvas->cd(5);
939   fMultiplicityESDCorrected->SetLineColor(kBlue);
940   fMultiplicityESDCorrected->Draw();
941   fMultiplicityMC->Draw("SAME");
942   fMultiplicityESD->Draw("SAME");
943
944   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
945   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
946   legend->AddEntry(fMultiplicityMC, "MC");
947   legend->AddEntry(fMultiplicityESD, "ESD");
948   legend->Draw();
949
950   canvas->cd(6);
951   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
952   fMultiplicityESDCorrectedRebinned->Draw();
953   fMultiplicityMC->Draw("SAME");
954
955   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
956   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
957   legend->AddEntry(fMultiplicityMC, "MC");
958   legend->Draw();
959
960   canvas->SaveAs("MultiplicityMC.gif");
961 }
962
963 void MultiplicityESD()
964 {
965   TFile* file = TFile::Open("multiplicityESD.root");
966
967   if (!file)
968   {
969     printf("multiplicityESD.root could not be opened.\n");
970     return;
971   }
972
973   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
974
975   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
976
977   fMultiplicityESD->Draw();
978 }
979
980 void drawPlots(Int_t max)
981 {
982   gMax = max;
983
984   ptCutoff();
985   TriggerBias();
986   VtxRecon();
987   Track2Particle2D();
988   Track2Particle3D();
989   ptSpectrum();
990   dNdEta();
991 }
992
993 void drawPlots()
994 {
995   drawPlots(5);
996   drawPlots(2);
997 }