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