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