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