adding flag for non-field data
[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("libVMC");
26   gSystem->Load("libTree");
27   gSystem->Load("libSTEERBase");
28   gSystem->Load("libESD");
29   gSystem->Load("libAOD");
30   gSystem->Load("libANALYSIS");
31   gSystem->Load("libANALYSISalice");
32   gSystem->Load("libPWG0base");
33 }
34
35 void SetRanges(TAxis* axis)
36 {
37   if (strcmp(axis->GetTitle(), "#eta") == 0)
38     axis->SetRangeUser(-1.4999, 1.4999);
39   if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0 || strcmp(axis->GetTitle(), "p_{T} (GeV/c)") == 0)
40   {
41     axis->SetRangeUser(0, 4.9999);
42     axis->SetTitle("p_{T} (GeV/c)");
43   }
44   if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0 || strcmp(axis->GetTitle(), "vtx z (cm)") == 0)
45   {
46     axis->SetRangeUser(-15, 14.9999);
47     axis->SetTitle("vtx-z (cm)");
48   }
49   if (strcmp(axis->GetTitle(), "Ntracks") == 0)
50     axis->SetRangeUser(0, 99.9999);
51 }
52
53 void SetRanges(TH1* hist)
54 {
55   SetRanges(hist->GetXaxis());
56   SetRanges(hist->GetYaxis());
57   SetRanges(hist->GetZaxis());
58 }
59
60
61 void Prepare3DPlot(TH3* hist)
62 {
63   hist->GetXaxis()->SetTitleOffset(1.5);
64   hist->GetYaxis()->SetTitleOffset(1.5);
65   hist->GetZaxis()->SetTitleOffset(1.5);
66
67   hist->SetStats(kFALSE);
68 }
69
70 void Prepare2DPlot(TH2* hist)
71 {
72   hist->SetStats(kFALSE);
73   hist->GetYaxis()->SetTitleOffset(1.4);
74
75   hist->SetMinimum(0);
76   hist->SetMaximum(gMax);
77
78   SetRanges(hist);
79 }
80
81 void Prepare1DPlot(TH1* hist)
82 {
83   hist->SetLineWidth(2);
84   hist->SetStats(kFALSE);
85
86   hist->GetXaxis()->SetLabelOffset(0.02);
87   hist->GetXaxis()->SetTitleOffset(1.3);
88   hist->GetYaxis()->SetTitleOffset(1.3);
89
90   SetRanges(hist);
91 }
92
93 void InitPad()
94 {
95   gPad->Range(0, 0, 1, 1);
96   gPad->SetLeftMargin(0.15);
97   //gPad->SetRightMargin(0.05);
98   //gPad->SetTopMargin(0.13);
99   gPad->SetBottomMargin(0.12);
100
101   gPad->SetGridx();
102   gPad->SetGridy();
103 }
104
105 void InitPadCOLZ()
106 {
107   gPad->Range(0, 0, 1, 1);
108   gPad->SetRightMargin(0.15);
109   gPad->SetLeftMargin(0.12);
110   gPad->SetTopMargin(0.05);
111
112   gPad->SetGridx();
113   gPad->SetGridy();
114 }
115
116 // --- end of helpers --- begin functions ---
117
118 void DrawOverview(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
119 {
120   loadlibs();
121   if (!TFile::Open(fileName))
122     return;
123
124   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
125   if (!dNdEtaCorrection->LoadHistograms())
126     return;
127
128   dNdEtaCorrection->Finish();
129
130   dNdEtaCorrection->DrawOverview();
131 }
132
133 void PrintInfo(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
134 {
135   loadlibs();
136   if (!TFile::Open(fileName))
137     return;
138
139   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
140   if (!dNdEtaCorrection->LoadHistograms())
141     return;
142
143   dNdEtaCorrection->Finish();
144
145   for (Int_t i=AlidNdEtaCorrection::kTrack2Particle; i<=AlidNdEtaCorrection::kND; i++)
146   {
147     Printf("Correction %d", i);
148     dNdEtaCorrection->GetCorrection(i)->PrintInfo(0);
149     return;
150   }
151 }
152
153 void PrintAllInfos()
154 {
155   PrintInfo();
156
157   Printf("RAW ESD");
158   TFile::Open("analysis_esd_raw.root");
159   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
160   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
161   fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
162
163   const Int_t num = 3;
164   const char* results[] = { "dndeta", "dndetaTr", "dndetaTrVtx" };
165
166   TFile::Open("analysis_esd.root");
167   for (Int_t i=0; i<num; i++)
168   {
169     Printf("CORRECTED %s", results[i]);
170     dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
171     fdNdEtaAnalysis->LoadHistograms(results[i]);
172     fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
173   }
174 }  
175
176 void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
177 {
178   gSystem->Load("libPWG0base");
179
180   TFile::Open(esdFile);
181   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
182   fdNdEtaAnalysisESD->LoadHistograms();
183
184   TFile::Open(mcFile);
185   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
186   fdNdEtaAnalysisMC->LoadHistograms();
187   //fdNdEtaAnalysisMC->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
188
189   for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
190     fdNdEtaAnalysisESD->GetdNdEtaPtCutOffCorrectedHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(i));
191
192   fdNdEtaAnalysisESD->DrawHistograms();
193 }
194
195 void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
196 {
197   gSystem->Load("libPWG0base");
198
199   const char* ESDfolder = 0;
200
201   if (plot == 0) // all
202     ESDfolder = "dndeta";
203   else if (plot == 1) // mb
204     ESDfolder = "dndeta_mb";
205   else if (plot == 2) // mb vtx
206     ESDfolder = "dndeta_mbvtx";
207
208   TFile::Open("analysis_esd.root");
209   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
210   fdNdEtaAnalysisESD->LoadHistograms();
211
212   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
213   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
214
215   TH2F* hist = 0;
216
217   if (plot == 0) // all
218     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
219   else if (plot == 1) // mb
220     hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
221   else if (plot == 2) // mb vtx
222     hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
223
224   TH1* proj = hist->ProjectionX();
225
226   TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
227   for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
228   {
229     Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
230     if (value != 0)
231     {
232       printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
233       vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
234     }
235   }
236
237   new TCanvas;
238   vertex->DrawCopy();
239 }
240
241 void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
242 {
243   loadlibs();
244
245   TFile::Open("analysis_esd.root");
246   dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
247   fdNdEtaAnalysisESD->LoadHistograms();
248
249   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
250   dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
251
252   //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
253   //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
254
255   TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
256   TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
257
258   TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
259
260   new TCanvas;
261   histESD->Draw();
262
263   new TCanvas;
264   histCorr->Draw();
265
266   for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
267     for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
268       for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
269       {
270         Float_t value1 = histESD->GetBinContent(x, y, z);
271         Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
272
273         if (value2 > 0 && value1 > 0)
274         {
275           printf("%f %f %f\n", value1, value2, value1 / value2);
276           diff->Fill(value1 / value2);
277         }
278       }
279
280   new TCanvas;
281   diff->Draw();
282 }
283
284 Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
285 {
286   Double_t avgMC = 0;
287   Double_t avgESD = 0;
288   for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
289   {
290     avgMC += histMC->GetBinContent(bin);
291     avgESD += histESD->GetBinContent(bin);
292   }
293   Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
294
295   avgMC /= nBins;
296   avgESD /= nBins;
297
298   // deviation when integrate in |eta| < 0.8 between mc and esd
299   Double_t diffFullRange = (avgMC - avgESD) / avgMC;
300
301   Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
302
303   return diffFullRange;
304 }
305
306 void dNdEtaNoResolution()
307 {
308   loadlibs();
309
310   TFile::Open("correction_map.root");
311
312   const char* correctionMapFolder = "dndeta_correction";
313   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
314   dNdEtaCorrection->LoadHistograms();
315   dNdEtaCorrection->GetTrack2ParticleCorrection()->PrintInfo(0.3);
316
317   TFile::Open("analysis_mc.root");
318   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
319   fdNdEtaAnalysis->LoadHistograms();
320   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, "ESD (no resolution effect) -> MB with vertex");
321   fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0)->SetMarkerStyle(21);
322
323   TFile::Open("analysis_mc.root");
324   dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
325   fdNdEtaAnalysisMC->LoadHistograms();
326   fdNdEtaAnalysisMC->Finish(0, 0, AlidNdEtaCorrection::kNone, "MC: MB with vertex");
327
328   DrawdNdEtaRatio(fdNdEtaAnalysis->GetdNdEtaPtCutOffCorrectedHistogram(0), fdNdEtaAnalysisMC->GetdNdEtaPtCutOffCorrectedHistogram(0), "MB with vertex (no resolution effect)", 3);
329 }
330
331 TH1* GetMCHist(const char* folder, Float_t ptCut, const char* tag)
332 {
333   loadlibs();
334   
335   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(folder, folder);
336   fdNdEtaAnalysis->LoadHistograms();
337   fdNdEtaAnalysis->Finish(0, ptCut, AlidNdEtaCorrection::kNone, tag);
338   return fdNdEtaAnalysis->GetdNdEtaHistogram(0);
339 }
340
341 void dNdEtaFinal(Bool_t spd = kTRUE)
342 {
343   TFile* file = TFile::Open("analysis_esd.root");
344   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
345   TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
346   Prepare1DPlot(histESD);
347   Prepare1DPlot(histESDnsd);
348   
349   TCanvas* canvas = new TCanvas("dNdEtaFinal", "dNdEtaFinal", 600, 600);
350   gPad->SetTopMargin(0.05);
351   gPad->SetRightMargin(0.05);
352   gPad->SetLeftMargin(0.12);
353   gPad->SetBottomMargin(0.12);
354   gPad->SetGridx();
355   gPad->SetGridy();
356   
357   Float_t etaMax = 1.9;
358   Float_t histMax = 1.39;
359   Float_t systErrorValue = 0.023;
360   Float_t systErrorNSDValue = 0.081;
361   if (!spd)
362   {
363     //etaMax = 1.5;
364     histMax = 0.99;
365     systErrorValue = 0.043;
366     systErrorNSDValue = 0.088;
367   }
368   
369   dummy = new TH2F("dummy", ";#eta;dN_{ch}/d#eta", 100, -etaMax, etaMax, 100, 3, 8);
370   dummy->SetStats(0);
371   dummy->GetYaxis()->SetTitleOffset(1.3);
372   
373   histESD->SetMarkerStyle(20);
374   histESDnsd->SetMarkerStyle(21);
375   histESDnsd->SetMarkerColor(4);
376   histESDnsd->SetLineColor(4);
377   histESD->SetMarkerSize(1.5);
378   histESDnsd->SetMarkerSize(1.5);
379   
380   histESD->GetXaxis()->SetRangeUser(-histMax, histMax);
381   histESDnsd->GetXaxis()->SetRangeUser(-histMax, histMax);
382   
383   legend = new TLegend(0.3, 0.2, 0.78, 0.4);
384   legend->SetFillColor(0);
385   legend->SetTextSize(0.04);
386   legend->AddEntry(histESD, "Inelastic events", "P");
387   legend->AddEntry(histESDnsd, "NSD events", "P");
388   
389   dummy->Draw();
390   
391   // syst errors.
392   TH1* systError = (TH1*) histESD->Clone("systError");
393   for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
394     systError->SetBinError(i, systError->GetBinContent(i) * systErrorValue);
395   // change error drawing style
396   systError->SetFillColor(15);    
397   systError->DrawCopy("SAME E2 ][");
398   
399   // syst error NSD
400   for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
401   {
402     systError->SetBinContent(i, histESDnsd->GetBinContent(i));
403     systError->SetBinError(i, systError->GetBinContent(i) * systErrorNSDValue);
404   }
405   systError->DrawCopy("SAME E2 ][");
406   
407   histESD->Draw("SAME");
408   histESDnsd->Draw("SAME");
409   legend->Draw();  
410   
411   canvas->SaveAs(Form("%s_dndeta_final.eps", (spd) ? "spd" : "tpc"));
412 }
413
414 void dNdEtaPythiaPhojet()
415 {
416   // evtl. deactivate acceptance maps in dNdEtaAnalysis.cxx
417
418   loadlibs();
419
420   TH1* hist[4];
421   
422   TFile::Open("LHC08c11_10TeV_0.5T/mb1/spd/analysis_mc.root");
423   hist[0] =         (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
424   hist[1] =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
425
426   TFile::Open("LHC08c15_10TeV_0.5T_Phojet/mb1/spd/analysis_mc.root");
427   hist[2] =         (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMCPhojet");
428   hist[3] =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsdPhojet");
429   
430   file = TFile::Open("pythia_phojet_dndeta.root", "RECREATE");
431   for (Int_t i=0; i<4; i++)
432     hist[i]->Write();
433   file->Close();
434 }
435  
436 void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
437 {
438   loadlibs();
439
440   TFile* file = TFile::Open("analysis_esd.root");
441   
442   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
443   fdNdEtaAnalysis->LoadHistograms("dndeta");
444   
445   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
446   TH1* histESDnsd = (TH1*) file->Get("dndetaNSD/dNdEta_corrected");
447   TH1* histESDnsdNoPt = (TH1*) file->Get("dndetaNSD/dNdEta");
448   TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
449   TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
450   TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
451   TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
452   TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
453   TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
454
455   Prepare1DPlot(histESD);
456   Prepare1DPlot(histESDnsd);
457   Prepare1DPlot(histESDMB);
458   Prepare1DPlot(histESDMBVtx);
459
460   Prepare1DPlot(histESDNoPt);
461   Prepare1DPlot(histESDMBNoPt);
462   Prepare1DPlot(histESDMBVtxNoPt);
463   Prepare1DPlot(histESDMBTracksNoPt);
464
465   histESD->SetLineWidth(0);
466   histESDnsd->SetLineWidth(0);
467   histESDMB->SetLineWidth(0);
468   histESDMBVtx->SetLineWidth(0);
469
470   histESDNoPt->SetLineWidth(0);
471   histESDMBNoPt->SetLineWidth(0);
472   histESDMBVtxNoPt->SetLineWidth(0);
473
474   histESD->SetMarkerColor(1);
475   histESDnsd->SetMarkerColor(6);
476   histESDMB->SetMarkerColor(2);
477   histESDMBVtx->SetMarkerColor(4);
478
479   histESD->SetLineColor(1);
480   histESDnsd->SetLineColor(6);
481   histESDMB->SetLineColor(2);
482   histESDMBVtx->SetLineColor(4);
483
484   histESDNoPt->SetMarkerColor(1);
485   histESDMBNoPt->SetMarkerColor(2);
486   histESDMBVtxNoPt->SetMarkerColor(4);
487   histESDMBTracksNoPt->SetMarkerColor(3);
488
489   histESD->SetMarkerStyle(20);
490   histESDnsd->SetMarkerStyle(29);
491   histESDMB->SetMarkerStyle(21);
492   histESDMBVtx->SetMarkerStyle(22);
493
494   histESDNoPt->SetMarkerStyle(20);
495   histESDMBNoPt->SetMarkerStyle(21);
496   histESDMBVtxNoPt->SetMarkerStyle(22);
497   histESDMBTracksNoPt->SetMarkerStyle(23);
498   
499   Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 0.89 : 1.79;
500   Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC) ? 1.2 : 2.3;
501   //Float_t etaLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 0.89 : 1.39;
502   //Float_t etaPlotLimit = (fdNdEtaAnalysis->GetAnalysisMode() == AliPWG0Helper::kTPC) ? 1.2 : 1.9;
503
504   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
505   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
506   histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
507   histESDnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
508
509   histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
510   histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
511   histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
512   histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
513
514   Float_t max = TMath::Max(histESDMBVtx->GetMaximum(), histESDMB->GetMaximum());
515   max = TMath::Max(max, histESD->GetMaximum());
516   max = TMath::Max(max, histESDMBTracksNoPt->GetMaximum());
517
518   TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.4);
519   legend->SetFillColor(0);
520   legend->AddEntry(histESDMBVtx, "Triggered, vertex");
521   legend->AddEntry(histESDMB, "Triggered");
522   legend->AddEntry(histESD, "All events");
523
524   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 2.1, max * 1.1);
525   Prepare1DPlot(dummy);
526   dummy->SetStats(kFALSE);
527   dummy->SetXTitle("#eta");
528   dummy->SetYTitle("dN_{ch}/d#eta");
529   dummy->GetYaxis()->SetTitleOffset(1);
530
531   TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
532
533   dummy->DrawCopy();
534   histESDMBVtx->Draw("SAME");
535   histESDMB->Draw("SAME");
536   histESD->Draw("SAME");
537   legend->Draw();
538
539   if (save)
540   {
541     canvas->SaveAs("dNdEta1.gif");
542     canvas->SaveAs("dNdEta1.eps");
543   }
544
545   if (onlyESD)
546     return;
547
548   loadlibs();
549
550   TFile* file2 = TFile::Open("analysis_mc.root");
551
552   TH1* histMCTrVtx =       (TH1*) GetMCHist("dndetaTrVtx", -1, "MC: MB with trigger")->Clone("histMCTrVtx");
553   TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
554   
555   TH1* histMC =            (TH1*) GetMCHist("dndeta", -1, "MC: full inelastic")->Clone("histMC");
556   TH1* histMCTr =          (TH1*) GetMCHist("dndetaTr", -1, "MC: minimum bias")->Clone("histMCTr");
557   TH1* histMCnsd =         (TH1*) GetMCHist("dndetaNSD", -1, "MC: NSD")->Clone("histMCnsd");
558
559   TH1* histMCPtCut =       (TH1*) GetMCHist("dndeta", 0.21, "MC: full inelastic, pt cut")->Clone("histMCPtCut");
560   TH1* histMCTrPtCut =     (TH1*) GetMCHist("dndetaTr", 0.21, "MC: minimum bias, pt cut")->Clone("histMCTrPtCut");
561   TH1* histMCTrVtxPtCut =  (TH1*) GetMCHist("dndetaTrVtx", 0.21, "MC: MB with trigger, pt cut")->Clone("histMCTrVtxPtCut");
562   TH1* histMCnsdNoPt =     (TH1*) GetMCHist("dndetaNSD", 0.21, "MC: NSD, put cut")->Clone("histMCnsdNoPt");
563   TH1* histMCTracksPtCut = (TH1*) GetMCHist("dndetaTracks", 0.21, "MC: Tracks w/o resolution effect, pt cut")->Clone("histMCTracksPtCut");
564
565   Prepare1DPlot(histMC);
566   Prepare1DPlot(histMCnsd);
567   Prepare1DPlot(histMCTr);
568   Prepare1DPlot(histMCTrVtx);
569
570   Prepare1DPlot(histMCPtCut);
571   Prepare1DPlot(histMCTrPtCut);
572   Prepare1DPlot(histMCTrVtxPtCut);
573   Prepare1DPlot(histMCTracksPtCut);
574
575   histMC->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
576   histMCnsd->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
577   histMCTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
578   histMCTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
579
580   histMCPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
581   histMCTrPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
582   histMCTrVtxPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
583   histMCTracksPtCut->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
584
585   histMC->SetLineColor(1);
586   histMCnsd->SetLineColor(6);
587   histMCTr->SetLineColor(2);
588   histMCTrVtx->SetLineColor(4);
589
590   histMCPtCut->SetLineColor(1);
591   histMCTrPtCut->SetLineColor(2);
592   histMCTrVtxPtCut->SetLineColor(4);
593   if (histMCTracksPtCut)
594     histMCTracksPtCut->SetLineColor(3);
595
596   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
597
598   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
599   dummy2->GetYaxis()->SetRangeUser(0, max * 1.1);
600
601   dummy2->DrawCopy();
602   histMC->Draw("SAME");
603   histMCnsd->Draw("SAME");
604   histMCTr->Draw("SAME");
605   histMCTrVtx->Draw("SAME");
606   histESD->Draw("SAME");
607   histESDnsd->Draw("SAME");
608   histESDMB->Draw("SAME");
609   histESDMBVtx->Draw("SAME");
610   histESDNoPt->Draw("SAME");
611   histESDMBNoPt->Draw("SAME");
612   histESDMBVtxNoPt->Draw("SAME");
613   histESDMBTracksNoPt->Draw("SAME");
614   histMCPtCut->Draw("SAME");
615   histMCTrPtCut->Draw("SAME");
616   histMCTrVtxPtCut->Draw("SAME");
617   if (histMCTracksPtCut)
618     histMCTracksPtCut->Draw("SAME");
619
620   if (save)
621   {
622     canvas2->SaveAs("dNdEta2.gif");
623     canvas2->SaveAs("dNdEta2.eps");
624   }
625
626   TH1* ratio = (TH1*) DrawdNdEtaRatio(histESD, histMC, "full_inelastic", etaPlotLimit)->Clone();
627   TH1* ratioTr = (TH1*) DrawdNdEtaRatio(histESDMB, histMCTr, "triggered", etaPlotLimit)->Clone();
628   TH1* ratioTrVtx = (TH1*) DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit)->Clone();
629   TH1* ratioTrVtxNoPt = (TH1*) DrawdNdEtaRatio(histESDMBVtxNoPt, histMCTrVtxPtCut, "triggered_vertex_nopt", etaPlotLimit)->Clone();
630   TH1* ratioNSD = (TH1*) DrawdNdEtaRatio(histESDnsd, histMCnsd, "NSD", etaPlotLimit)->Clone();
631
632   // draw ratios of single steps
633   c7 = new TCanvas("all_ratios", "all_ratios", 600, 600);
634   c7->SetRightMargin(0.05);
635   c7->SetTopMargin(0.05);
636   c7->SetGridx();
637   c7->SetGridy();
638   
639   ratioTrVtxNoPt->SetMarkerStyle(20);
640   ratioTrVtx->SetMarkerStyle(21);
641   ratioTr->SetMarkerStyle(23);
642   ratio->SetMarkerStyle(22);
643   ratioNSD->SetMarkerStyle(26);
644   
645   ratioTrVtxNoPt->SetMarkerSize(2);
646   ratioTrVtx->SetMarkerSize(2);
647   ratioTr->SetMarkerSize(2);
648   ratio->SetMarkerSize(2);
649   ratioNSD->SetMarkerSize(2);
650   
651   ratioTrVtxNoPt->SetMarkerColor(1);
652   ratioTrVtx->SetMarkerColor(2);
653   ratioTr->SetMarkerColor(4);
654   ratio->SetMarkerColor(2);
655   ratioNSD->SetMarkerColor(1);
656   
657   ratioTrVtxNoPt->SetLineColor(1);
658   ratioTrVtx->SetLineColor(2);
659   ratioTr->SetLineColor(4);
660   ratio->SetLineColor(2);
661   ratioNSD->SetLineColor(1);
662   
663   legend7 = new TLegend(0.13, 0.7, 0.94, 0.9);
664   legend7->SetFillColor(0);
665   legend7->SetTextSize(0.035);
666   legend7->SetNColumns(2);
667   
668   flat = new TF1("flat", "-1", -5, 5);
669   ratioTrVtxNoPt->Add(flat);
670   ratioTrVtx->Add(flat);
671   ratioTr->Add(flat);
672   ratio->Add(flat);
673   ratioNSD->Add(flat);
674   
675   ratioTrVtxNoPt->Scale(100);
676   ratioTrVtx->Scale(100);
677   ratioTr->Scale(100);
678   ratio->Scale(100);
679   ratioNSD->Scale(100);
680   
681   ratio->Add(ratioTr, -1);
682   ratioNSD->Add(ratioTr, -1);
683   ratioTr->Add(ratioTrVtx, -1);
684   ratioTrVtx->Add(ratioTrVtxNoPt, -1);
685   
686   legend7->AddEntry(ratioTrVtxNoPt, "Track-to-particle", "P");
687   legend7->AddEntry(ratio, "Trigger-bias INEL", "P");
688   legend7->AddEntry(ratioTr, "Vertex-reconstruction", "P");
689   legend7->AddEntry(ratioNSD, "Trigger-bias NSD", "P");
690   if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
691     legend7->AddEntry(ratioTrVtx, "p_{T} cut-off", "P");
692   
693   TH1* dummy7 = new TH2F("dummy7", ";#eta;Deviation in %", 100, -etaPlotLimit, etaPlotLimit, 100, -5, 7);
694   dummy7->SetStats(0);
695   dummy7->Draw();
696   
697   ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
698   ratioTr->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
699   ratioTrVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
700   ratioTrVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
701   ratioNSD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
702   
703   ratio->Draw("HIST EP SAME");
704   ratioTr->Draw("HIST EP SAME");
705   if ((fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kFieldOn) && (fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPC || fdNdEtaAnalysis->GetAnalysisMode() & AliPWG0Helper::kTPCITS))
706     ratioTrVtx->Draw("HIST EP SAME");
707   ratioTrVtxNoPt->Draw("HIST EP SAME");
708   ratioNSD->Draw("HIST EP SAME");
709   legend7->Draw();
710   
711   c7->SaveAs("ratios.eps");
712
713   new TCanvas;
714   dummy2->DrawCopy();
715   histMCnsd->Draw("SAME");
716   histESDnsd->Draw("SAME");
717
718   ratio = (TH1*) histMC->Clone("ratio");
719   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
720
721   ratio->Divide(histESD);
722   ratioNoPt->Divide(histESDNoPt);
723
724   ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
725
726   ratio->SetLineColor(1);
727   ratioNoPt->SetLineColor(2);
728
729   Double_t average = 0;       // average deviation from 1 in ratio (depends on the number of bins if statistical)
730   for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
731     average += TMath::Abs(ratio->GetBinContent(bin) - 1);
732   Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
733   average /= nBins;
734   Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
735
736   PrintIntegratedDeviation(histMC, histESD, "all events");
737   PrintIntegratedDeviation(histMCnsd, histESDnsd, "all events (NSD)");
738   PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
739   PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
740   PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
741   PrintIntegratedDeviation(histMCnsdNoPt, histESDnsdNoPt, "all events (NSD) (no pt corr)");
742   PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
743   PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
744   PrintIntegratedDeviation(histESD, histESDNoPt, "pt cut off correction");
745
746   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 600, 600);
747   canvas3->Range(0, 0, 1, 1);
748   //canvas3->Divide(1, 2, 0, 0);
749
750   //canvas3->cd(1);
751   TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.5, 0.98, 0.98);
752   pad1->SetTopMargin(0.05);
753   pad1->SetLeftMargin(0.13);
754   pad1->Draw();
755
756   TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.5);
757   pad2->SetLeftMargin(0.13);
758   pad2->Draw();
759
760   pad1->SetRightMargin(0.01);
761   pad2->SetRightMargin(0.01);
762
763   // no border between them
764   pad1->SetBottomMargin(0);
765   pad2->SetTopMargin(0);
766
767   pad1->cd();
768   pad1->SetGridx();
769   pad1->SetGridy();
770
771   legend->AddEntry(histMC, "MC prediction");
772
773   dummy->GetXaxis()->SetLabelSize(0.08);
774   dummy->GetYaxis()->SetLabelSize(0.08);
775   dummy->GetXaxis()->SetTitleSize(0.08);
776   dummy->GetYaxis()->SetTitleSize(0.08);
777   dummy->GetYaxis()->SetTitleOffset(0.8);
778   dummy->DrawCopy();
779   histESDMBVtx->Draw("SAME");
780   histESDMB->Draw("SAME");
781   histESD->Draw("SAME");
782   histMC->Draw("SAME");
783
784   legend->SetTextSize(0.08);
785   legend->Draw();
786
787   pad2->cd();
788   pad2->SetBottomMargin(0.15);
789   //pad2->SetGridx();
790   //pad2->SetGridy();
791
792   Float_t minR = 0.91; //TMath::Min(0.961, ratio->GetMinimum() * 0.95);
793   Float_t maxR = 1.09; //TMath::Max(1.049, ratio->GetMaximum() * 1.05);
794
795   TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
796   dummy3.SetStats(kFALSE);
797   for (Int_t i=1; i<=100; ++i)
798     dummy3.SetBinContent(i, 1);
799   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
800   dummy3.SetLineWidth(2);
801   dummy3.GetXaxis()->SetLabelSize(0.08);
802   dummy3.GetYaxis()->SetLabelSize(0.08);
803   dummy3.GetXaxis()->SetTitleSize(0.08);
804   dummy3.GetYaxis()->SetTitleSize(0.08);
805   dummy3.GetYaxis()->SetTitleOffset(0.8);
806   dummy3.DrawCopy();
807
808   ratio->Draw("SAME");
809
810   //pad2->Draw();
811
812   canvas3->Modified();
813
814   if (save)
815   {
816     canvas3->SaveAs("dNdEta.gif");
817     canvas3->SaveAs("dNdEta.eps");
818   }
819
820   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
821
822   ratio->Draw();
823   ratioNoPt->Draw("SAME");
824
825   TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.9);
826   legend->SetFillColor(0);
827   legend->AddEntry(ratio, "mc/esd");
828   legend->AddEntry(ratioNoPt, "mc/esd, not pt cut off corrected");
829   legend->Draw();
830 }
831
832 TH1* DrawdNdEtaRatio(TH1* corr, TH1* mc, const char* name, Float_t etaPlotLimit)
833 {
834   TCanvas* canvas3 = new TCanvas(name, name, 600, 600);
835   canvas3->Range(0, 0, 1, 1);
836
837   TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
838   pad1->Draw();
839
840   TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
841   pad2->Draw();
842
843   pad1->SetRightMargin(0.01);
844   pad2->SetRightMargin(0.01);
845   pad1->SetTopMargin(0.05);
846   pad1->SetLeftMargin(0.13);
847   pad2->SetLeftMargin(0.13);
848   pad2->SetBottomMargin(0.15);
849   
850   // no border between them
851   pad1->SetBottomMargin(0);
852   pad2->SetTopMargin(0);
853
854   pad1->cd();
855   pad1->SetGridx();
856   pad1->SetGridy();
857
858   TLegend* legend = new TLegend(0.35, 0.05, 0.75, 0.3);
859   legend->SetFillColor(0);
860   legend->AddEntry(corr, "Corrected");
861   legend->AddEntry(mc, "MC prediction");
862   legend->SetTextSize(0.08);
863
864   TH2F* dummy = new TH2F("dummy", "", 100, -etaPlotLimit, etaPlotLimit, 1000, 2.7, corr->GetMaximum() * 1.1);
865   Prepare1DPlot(dummy);
866   dummy->SetStats(kFALSE);
867   dummy->SetXTitle("#eta");
868   dummy->SetYTitle("dN_{ch}/d#eta");
869   dummy->GetYaxis()->SetTitleOffset(1);
870
871   dummy->GetXaxis()->SetLabelSize(0.08);
872   dummy->GetYaxis()->SetLabelSize(0.08);
873   dummy->GetXaxis()->SetTitleSize(0.08);
874   dummy->GetYaxis()->SetTitleSize(0.08);
875   dummy->GetYaxis()->SetTitleOffset(0.8);
876   dummy->DrawCopy();
877
878   corr->Draw("SAME");
879   mc->Draw("SAME");
880
881   legend->Draw();
882
883   pad2->cd();
884   pad2->SetBottomMargin(0.15);
885   //pad2->SetGridx();
886   //pad2->SetGridy();
887
888   TH1* ratio = (TH1*) mc->Clone("ratio");
889   ratio->Divide(corr);
890
891   Float_t minR = TMath::Min(0.91, ratio->GetMinimum() * 0.95);
892   Float_t maxR = TMath::Max(1.09, ratio->GetMaximum() * 1.05);
893
894   TH1F dummy3("dummy3", ";#eta;Ratio: MC / corr", 100, -etaPlotLimit, etaPlotLimit);
895   dummy3.SetStats(kFALSE);
896   for (Int_t i=1; i<=100; ++i)
897         dummy3.SetBinContent(i, 1);
898   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
899   dummy3.SetLineWidth(2);
900   dummy3.GetXaxis()->SetLabelSize(0.08);
901   dummy3.GetYaxis()->SetLabelSize(0.08);
902   dummy3.GetXaxis()->SetTitleSize(0.08);
903   dummy3.GetYaxis()->SetTitleSize(0.08);
904   dummy3.GetYaxis()->SetTitleOffset(0.8);
905   dummy3.DrawCopy();
906
907   ratio->Draw("SAME");
908
909   canvas3->Modified();
910
911   return ratio;
912 }
913
914 void ptSpectrum()
915 {
916   TFile* file = TFile::Open("analysis_esd.root");
917   TH1* histESD = (TH1*) file->Get("dndeta/dndeta_pt");
918
919   TFile* file2 = TFile::Open("analysis_mc.root");
920   TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_pt");
921
922   TCanvas* canvas = new TCanvas("ptSpectrum", "ptSpectrum", 500, 500);
923   InitPad();
924   gPad->SetLogy();
925
926   Prepare1DPlot(histMC);
927   Prepare1DPlot(histESD);
928
929   histESD->SetTitle("");
930   histESD->GetXaxis()->SetTitle("p_{T} [GeV/c]");
931   histESD->GetYaxis()->SetTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
932
933   histMC->SetLineColor(kBlue);
934   histESD->SetLineColor(kRed);
935
936   histESD->GetYaxis()->SetTitleOffset(1.5);
937   histESD->GetXaxis()->SetRangeUser(0, 4.9999);
938
939   histESD->SetMaximum(TMath::Max(histESD->GetMaximum(), histMC->GetMaximum()) * 2);
940
941   histESD->Draw();
942   histMC->Draw("SAME");
943
944   canvas->SaveAs("ptSpectrum.gif");
945   canvas->SaveAs("ptSpectrum.eps");
946 }
947
948 void TriggerBiasVtxRecon(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
949 {
950   gSystem->Load("libPWG0base");
951
952   TFile::Open(fileName);
953   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
954   dNdEtaCorrection->LoadHistograms();
955
956   TH2* corrTrigger = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram();
957   TH2* corrVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram();
958
959   Prepare2DPlot(corrTrigger);
960   corrTrigger->SetTitle("b) Trigger bias correction");
961
962   Prepare2DPlot(corrVtx);
963   corrVtx->SetTitle("a) Vertex reconstruction correction");
964
965   corrTrigger->GetYaxis()->SetTitle("Multiplicity");
966   corrVtx->GetYaxis()->SetTitle("Multiplicity");
967
968   TCanvas* canvas = new TCanvas("TriggerBiasVtxRecon", "TriggerBiasVtxRecon", 1000, 500);
969   canvas->Divide(2, 1);
970
971   canvas->cd(1);
972   InitPadCOLZ();
973   corrVtx->DrawCopy("COLZ");
974
975   canvas->cd(2);
976   InitPadCOLZ();
977   corrTrigger->DrawCopy("COLZ");
978
979   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.gif", gMax));
980   canvas->SaveAs(Form("TriggerBiasVtxRecon_%d.eps", gMax));
981
982   canvas = new TCanvas("TriggerBiasVtxReconZoom", "TriggerBiasVtxReconZoom", 1000, 500);
983   canvas->Divide(2, 1);
984
985   corrTrigger->GetYaxis()->SetRangeUser(0, 5);
986   corrVtx->GetYaxis()->SetRangeUser(0, 5);
987
988   canvas->cd(1);
989   InitPadCOLZ();
990   corrVtx->DrawCopy("COLZ");
991
992   canvas->cd(2);
993   InitPadCOLZ();
994   corrTrigger->DrawCopy("COLZ");
995
996   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.gif", gMax));
997   canvas->SaveAs(Form("TriggerBiasVtxReconZoom_%d.eps", gMax));
998 }
999
1000 void TriggerBias(const char* fileName = "correction_map.root")
1001 {
1002   TFile* file = TFile::Open(fileName);
1003
1004   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_trigger"));
1005
1006   Prepare2DPlot(corr);
1007   corr->SetTitle("Trigger bias correction");
1008
1009   TCanvas* canvas = new TCanvas("TriggerBias", "TriggerBias", 500, 500);
1010   InitPadCOLZ();
1011   corr->DrawCopy("COLZ");
1012
1013   canvas->SaveAs(Form("TriggerBias_%d.gif", gMax));
1014   canvas->SaveAs(Form("TriggerBias_%d.eps", gMax));
1015
1016   corr->GetYaxis()->SetRangeUser(0, 5);
1017
1018   canvas = new TCanvas("TriggerBiasZoom", "TriggerBiasZoom", 500, 500);
1019   InitPadCOLZ();
1020   corr->DrawCopy("COLZ");
1021
1022   canvas->SaveAs(Form("TriggerBiasZoom_%d.gif", gMax));
1023   canvas->SaveAs(Form("TriggerBiasZoom_%d.eps", gMax));
1024 }
1025
1026 void TriggerBias1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1027 {
1028   gSystem->Load("libPWG0base");
1029
1030   TFile* file = TFile::Open(fileName);
1031   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1032   dNdEtaCorrection->LoadHistograms();
1033
1034   TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("x");
1035   TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1036
1037   TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 1000, 500);
1038   canvas->Divide(2, 1);
1039
1040   canvas->cd(1);
1041   InitPad();
1042
1043   Prepare1DPlot(hist);
1044   hist->SetTitle("");
1045   hist->GetYaxis()->SetTitle("correction factor");
1046   hist->GetYaxis()->SetRangeUser(1, 1.5);
1047   hist->GetYaxis()->SetTitleOffset(1.6);
1048   hist->Draw();
1049
1050   canvas->cd(2);
1051   InitPad();
1052
1053   Prepare1DPlot(hist2);
1054   hist2->SetTitle("");
1055   hist2->GetYaxis()->SetTitle("correction factor");
1056   hist2->GetXaxis()->SetRangeUser(0, 5);
1057   hist2->GetYaxis()->SetTitleOffset(1.6);
1058   hist2->GetXaxis()->SetTitle("multiplicity");
1059   hist2->Draw();
1060
1061   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1062   pave->SetFillColor(0);
1063   pave->AddText("|z| < 10 cm");
1064   pave->Draw();
1065
1066   canvas->SaveAs("TriggerBias1D.eps");
1067 }
1068
1069 void VtxRecon()
1070 {
1071   TFile* file = TFile::Open("correction_map.root");
1072
1073   TH2* corr = dynamic_cast<TH2*> (file->Get("dndeta_correction/corr_dndeta_correction_vtxReco"));
1074
1075   Prepare2DPlot(corr);
1076   corr->SetTitle("Vertex reconstruction correction");
1077
1078   TCanvas* canvas = new TCanvas("VtxRecon", "VtxRecon", 500, 500);
1079   InitPadCOLZ();
1080   corr->DrawCopy("COLZ");
1081
1082   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
1083   canvas->SaveAs(Form("VtxRecon_%d.eps", gMax));
1084
1085   corr->GetYaxis()->SetRangeUser(0, 5);
1086
1087   canvas = new TCanvas("VtxReconZoom", "VtxReconZoom", 500, 500);
1088   InitPadCOLZ();
1089   corr->DrawCopy("COLZ");
1090
1091   canvas->SaveAs(Form("VtxReconZoom_%d.gif", gMax));
1092   canvas->SaveAs(Form("VtxReconZoom_%d.eps", gMax));
1093 }
1094
1095 void VtxRecon1D(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction")
1096 {
1097   gSystem->Load("libPWG0base");
1098
1099   TFile* file = TFile::Open(fileName);
1100   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1101   dNdEtaCorrection->LoadHistograms();
1102
1103   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("x");
1104   TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1105
1106   TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
1107   canvas->Divide(2, 1);
1108
1109   canvas->cd(1);
1110   InitPad();
1111
1112   Prepare1DPlot(hist);
1113   hist->SetTitle("");
1114   hist->GetYaxis()->SetTitle("correction factor");
1115   hist->GetYaxis()->SetRangeUser(1, 1.8);
1116   hist->GetYaxis()->SetTitleOffset(1.6);
1117   hist->DrawCopy();
1118
1119   canvas->cd(2);
1120   InitPad();
1121
1122   Prepare1DPlot(hist2);
1123   hist2->SetTitle("");
1124   hist2->GetYaxis()->SetTitle("correction factor");
1125   hist2->GetXaxis()->SetRangeUser(0, 20);
1126   hist2->GetYaxis()->SetTitleOffset(1.6);
1127   hist2->GetXaxis()->SetTitle("multiplicity");
1128   hist2->Draw();
1129
1130   TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
1131   pave->SetFillColor(0);
1132   pave->AddText("|z| < 10 cm");
1133   pave->Draw();
1134
1135   canvas->SaveAs("VtxRecon1D.eps");
1136
1137   Correction1DCreatePlots(fileName, folderName, 9.9, 2);
1138
1139   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1140   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1141
1142   Prepare1DPlot(corrX);
1143   Prepare1DPlot(corrZ);
1144
1145   corrX->GetYaxis()->SetTitleOffset(1.5);
1146   corrZ->GetYaxis()->SetTitleOffset(1.5);
1147
1148   corrX->SetTitle("a) z projection");
1149   corrZ->SetTitle("b) p_{T} projection");
1150
1151   corrX->GetYaxis()->SetTitle("Correction factor");
1152   corrZ->GetYaxis()->SetTitle("Correction factor");
1153
1154   corrZ->GetXaxis()->SetRangeUser(0.11, 9.9);
1155
1156   TString canvasName;
1157   canvasName.Form("VtxRecon1D_Track");
1158   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
1159   canvas->Divide(2, 1);
1160
1161   canvas->cd(1);
1162   InitPad();
1163   corrX->DrawCopy();
1164
1165   canvas->cd(2);
1166   InitPad();
1167   gPad->SetLogx();
1168   corrZ->Draw();
1169
1170   canvas->SaveAs("VtxRecon1D_Track.eps");
1171   canvas->SaveAs("VtxRecon1D_Track.gif");
1172 }
1173
1174 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
1175 {
1176   gSystem->Load("libPWG0base");
1177
1178   TFile::Open(fileName);
1179   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1180   dNdEtaCorrection->LoadHistograms(fileName, "dndeta_correction");
1181
1182   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1183   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1184
1185   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1186   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1187
1188   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1189   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1190   gene->GetXaxis()->SetRangeUser(-10, 10);
1191   meas->GetXaxis()->SetRangeUser(-10, 10);
1192
1193   Float_t eff1 = gene->Integral() / meas->Integral();
1194   Float_t error1 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1195
1196   printf("Correction without pT cut: %f +- %f\n", eff1, error1);
1197
1198   gene->GetZaxis()->SetRangeUser(0.3, 10);
1199   meas->GetZaxis()->SetRangeUser(0.3, 10);
1200
1201   Float_t eff2 = gene->Integral() / meas->Integral();
1202   Float_t error2 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1203
1204   printf("Correction with pT cut: %f +- %f\n", eff2, error2);
1205
1206   gene->GetZaxis()->SetRangeUser(0.3, 1);
1207   meas->GetZaxis()->SetRangeUser(0.3, 1);
1208
1209   Float_t eff3 = gene->Integral() / meas->Integral();
1210   Float_t error3 = TMath::Sqrt(gene->Integral()) / meas->Integral();
1211
1212   printf("Correction with 0.3 < pT < 0.5: %f +- %f\n", eff3, error3);
1213 }
1214
1215 void Correction1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType = 0, Int_t correctionType2 = -1)
1216 {
1217   if (correctionType2 == -1)
1218     correctionType2 = correctionType;
1219
1220   TFile::Open(fileName);
1221   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
1222   dNdEtaCorrection->LoadHistograms();
1223
1224   TH3F* gene = dNdEtaCorrection->GetCorrection(correctionType)->GetTrackCorrection()->GetGeneratedHistogram();
1225   TH3F* meas = dNdEtaCorrection->GetCorrection(correctionType2)->GetTrackCorrection()->GetMeasuredHistogram();
1226
1227   gene->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1228   meas->GetZaxis()->SetRangeUser(0.3, upperPtLimit);
1229   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1230   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1231   AliPWG0Helper::CreateDividedProjections(gene, meas, "x", kFALSE);
1232   gene->GetYaxis()->SetRange(0, 0);
1233   meas->GetYaxis()->SetRange(0, 0);
1234
1235   gene->GetXaxis()->SetRangeUser(-9.9, 9.9);
1236   meas->GetXaxis()->SetRangeUser(-9.9, 9.9);
1237   AliPWG0Helper::CreateDividedProjections(gene, meas, "y", kFALSE);
1238   gene->GetZaxis()->SetRange(0, 0);
1239   meas->GetZaxis()->SetRange(0, 0);
1240
1241   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1242   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1243   AliPWG0Helper::CreateDividedProjections(gene, meas, "z", kFALSE);
1244 }
1245
1246 TCanvas* Correction1D(Int_t correctionType = 0, const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9, Int_t correctionType2 = -1)
1247 {
1248   gSystem->Load("libPWG0base");
1249
1250   Correction1DCreatePlots(fileName, folder, upperPtLimit, correctionType, correctionType2);
1251
1252   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1253   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1254   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1255
1256   Prepare1DPlot(corrX);
1257   Prepare1DPlot(corrY);
1258   Prepare1DPlot(corrZ);
1259
1260   /*
1261   corrX->SetTitle("a) z projection");
1262   corrY->SetTitle("b) #eta projection");
1263   corrZ->SetTitle("c) p_{T} projection");
1264   */
1265   
1266   corrX->SetTitle("");
1267   corrY->SetTitle("");
1268   corrZ->SetTitle("");
1269
1270   corrX->SetTitleSize(0.06, "xyz");
1271   corrX->SetLabelSize(0.06, "xyz");
1272   corrY->SetTitleSize(0.06, "xyz");
1273   corrY->SetLabelSize(0.06, "xyz");
1274   corrZ->SetTitleSize(0.06, "xyz");
1275   corrZ->SetLabelSize(0.06, "xyz");
1276
1277   corrX->GetYaxis()->SetTitle("Correction factor");
1278   corrY->GetYaxis()->SetTitle("Correction factor");
1279   corrZ->GetYaxis()->SetTitle("Correction factor");
1280   //corrX->GetYaxis()->SetTitleOffset(1.7);
1281   //corrY->GetYaxis()->SetTitleOffset(1.7);
1282   //corrZ->GetYaxis()->SetTitleOffset(1.7);
1283   corrX->GetYaxis()->SetRangeUser(0.8, 1.5);
1284   corrY->GetYaxis()->SetRangeUser(0.8, 1.5);
1285   corrZ->GetYaxis()->SetRangeUser(0.8, 1.5);
1286
1287   corrZ->GetXaxis()->SetRangeUser(0.11, upperPtLimit);
1288
1289   TString canvasName;
1290   canvasName.Form(Form("Correction1D_%d_%s_%f", correctionType, fileName, upperPtLimit));
1291   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1292   canvas->Divide(3, 1);
1293
1294   TLatex* Tl = new TLatex;
1295   Tl->SetTextSize(0.06);
1296   Tl->SetBit(TLatex::kTextNDC);
1297
1298   canvas->cd(1);
1299   InitPad();
1300   gPad->SetTopMargin(0.05);
1301   gPad->SetBottomMargin(0.15);
1302   corrX->DrawCopy();
1303   Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1304   Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
1305
1306   canvas->cd(2);
1307   InitPad();
1308   gPad->SetTopMargin(0.05);
1309   gPad->SetBottomMargin(0.15);
1310   corrY->Draw();
1311   Tl->DrawLatex(0.5, 0.88, "0.3 < p_{T} < 10");
1312   Tl->DrawLatex(0.5, 0.8, "|vtx-z| < 10 cm");
1313
1314   canvas->cd(3);
1315   InitPad();
1316   gPad->SetTopMargin(0.05);
1317   gPad->SetBottomMargin(0.15);
1318   gPad->SetLogx();
1319   corrZ->Draw();
1320   corrZ->GetXaxis()->SetLabelOffset(0.005);
1321   corrZ->GetXaxis()->SetTitleOffset(1.2);
1322   Tl->DrawLatex(0.5, 0.88, "|vtx-z| < 10 cm");
1323   Tl->DrawLatex(0.5, 0.8, "|#eta| < 0.8");
1324
1325   return canvas;
1326 }
1327
1328 void Track2Particle1D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction", Float_t upperPtLimit = 9.9)
1329 {
1330   gSystem->Load("libPWG0base");
1331
1332   Correction1DCreatePlots(fileName, folder, upperPtLimit, AlidNdEtaCorrection::kTrack2Particle);
1333
1334   TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_x_div_measured_x", folder, folder)));
1335   TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_y_div_measured_y", folder, folder)));
1336   TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("generated_z_div_measured_z", folder, folder)));
1337
1338   Prepare1DPlot(corrX);
1339   Prepare1DPlot(corrY);
1340   Prepare1DPlot(corrZ);
1341
1342   corrX->SetTitle("a) z projection");
1343   corrY->SetTitle("a) #eta projection");
1344   corrZ->SetTitle("b) p_{T} projection");
1345
1346   corrY->GetYaxis()->SetTitle("correction factor");
1347   corrZ->GetYaxis()->SetTitle("correction factor");
1348
1349   corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1350
1351   TString canvasName;
1352   canvasName.Form("Track2Particle1D_%s", folder);
1353   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1354   canvas->Divide(3, 1);
1355
1356   canvas->cd(1);
1357   InitPad();
1358   corrX->DrawCopy();
1359
1360   canvas->cd(2);
1361   InitPad();
1362   corrY->Draw();
1363
1364   canvas->cd(3);
1365   InitPad();
1366   corrZ->Draw();
1367
1368   canvas->SaveAs(Form("Track2Particle1D_%s_%f.gif", fileName, upperPtLimit));
1369   canvas->SaveAs(Form("Track2Particle1D_%s_%f.eps", fileName, upperPtLimit));
1370
1371   //TPaveText* pave = new TPaveText(-0.4, 1.35, 0.4, 1.45);
1372
1373   canvasName.Form("Track2Particle1D_%s_etapt", folder);
1374   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
1375   canvas->Divide(2, 1);
1376
1377   canvas->cd(1);
1378   InitPad();
1379   corrY->GetXaxis()->SetRangeUser(-0.99, 0.99);
1380   corrY->GetYaxis()->SetRangeUser(1, 1.5);
1381   corrY->GetYaxis()->SetTitleOffset(1.5);
1382   corrY->DrawCopy();
1383   TPaveText* pave = new TPaveText(0.3, 0.7, 0.7, 0.8, "NDC");
1384   pave->AddText("|z| < 10 cm");
1385   pave->AddText("0.3 GeV/c < p_{T} < 10 GeV/c");
1386   pave->Draw();
1387
1388   canvas->cd(2);
1389   InitPad();
1390   gPad->SetLogx();
1391   corrZ->GetYaxis()->SetRangeUser(1, 2.5);
1392   corrZ->GetXaxis()->SetRangeUser(0.101, upperPtLimit);
1393   corrZ->GetYaxis()->SetTitleOffset(1.5);
1394   corrZ->DrawCopy();
1395   pave = new TPaveText(0.5, 0.7, 0.8, 0.8, "NDC");
1396   pave->AddText("|z| < 10 cm");
1397   pave->AddText("|#eta| < 0.8");
1398   pave->Draw();
1399
1400   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.eps", fileName, upperPtLimit));
1401   canvas->SaveAs(Form("Track2Particle1D_etapt_%s_%f.gif", fileName, upperPtLimit));
1402 }
1403
1404 /*
1405 void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
1406 {
1407   gSystem->Load("libPWG0base");
1408
1409   // particle type
1410   for (Int_t particle=0; particle<4; ++particle)
1411   {
1412     TString dirName;
1413     dirName.Form("correction_%d", particle);
1414     Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
1415
1416     TString tmpx, tmpy, tmpz;
1417     tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
1418     tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
1419     tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
1420
1421     TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
1422     TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
1423     TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
1424
1425     Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
1426
1427     TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
1428     TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
1429     TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
1430
1431     posX->Divide(negX);
1432     posY->Divide(negY);
1433     posZ->Divide(negZ);
1434
1435     Prepare1DPlot(posX);
1436     Prepare1DPlot(posY);
1437     Prepare1DPlot(posZ);
1438
1439     Float_t min = 0.8;
1440     Float_t max = 1.2;
1441
1442     posX->SetMinimum(min);
1443     posX->SetMaximum(max);
1444     posY->SetMinimum(min);
1445     posY->SetMaximum(max);
1446     posZ->SetMinimum(min);
1447     posZ->SetMaximum(max);
1448
1449     posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
1450
1451     posX->GetYaxis()->SetTitleOffset(1.7);
1452     posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
1453     posY->GetYaxis()->SetTitleOffset(1.7);
1454     posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
1455     posZ->GetYaxis()->SetTitleOffset(1.7);
1456     posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
1457
1458     posZ->GetXaxis()->SetRangeUser(0, 1);
1459
1460     TString canvasName;
1461     canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
1462
1463     TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
1464     canvas->Divide(3, 1);
1465
1466     canvas->cd(1);
1467     InitPad();
1468     posX->DrawCopy();
1469
1470     canvas->cd(2);
1471     InitPad();
1472     posY->DrawCopy();
1473
1474     canvas->cd(3);
1475     InitPad();
1476     posZ->DrawCopy();
1477
1478     canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1479     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
1480   }
1481 }
1482 */
1483
1484 void CompareTrack2Particle1D(const char* file1, const char* file2, Float_t upperPtLimit = 9.9)
1485 {
1486   loadlibs();
1487
1488   const char* folderName = "dndeta_correction";
1489
1490   c = new TCanvas("CompareTrack2Particle1D", "CompareTrack2Particle1D", 1200, 400);
1491   c->Divide(3, 1);
1492
1493   for (Int_t fileId = 0; fileId < 2; fileId++)
1494   {
1495     const char* file = ((fileId == 0) ? file1 : file2);
1496     Correction1DCreatePlots(file, folderName, upperPtLimit, 1);
1497
1498     TH1* corr[3];
1499     corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
1500     corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y_div_measured_y"));
1501     corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
1502     /*corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x"))->Clone(Form("hist_x_%d", fileId));
1503     corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y"))->Clone(Form("hist_y_%d", fileId));
1504     corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z"))->Clone(Form("hist_z_%d", fileId));*/
1505
1506     for (Int_t i=0; i<3; i++)
1507     {
1508       c->cd(i+1);
1509       InitPad();
1510       corr[i]->GetYaxis()->SetRangeUser(0.8, 2);
1511       corr[i]->SetLineColor(fileId+1);
1512       corr[i]->DrawCopy((fileId == 0) ? "" : "SAME");
1513     }
1514   }
1515
1516   return;
1517
1518   c->SaveAs(Form("%s.gif", canvas->GetName()));
1519   c->SaveAs(Form("%s.eps", canvas->GetName()));
1520 }
1521
1522 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
1523 {
1524   TFile::Open(fileName);
1525   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
1526   dNdEtaCorrection->LoadHistograms();
1527
1528   TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
1529   TH3F* meas = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
1530
1531   gene->GetZaxis()->SetRangeUser(0.2, 10);
1532   meas->GetZaxis()->SetRangeUser(0.2, 10);
1533   AliPWG0Helper::CreateDividedProjections(gene, meas, "yx");
1534   gene->GetZaxis()->SetRange(0, 0);
1535   meas->GetZaxis()->SetRange(0, 0);
1536
1537   gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
1538   meas->GetYaxis()->SetRangeUser(-0.8, 0.8);
1539   AliPWG0Helper::CreateDividedProjections(gene, meas, "zx");
1540   gene->GetYaxis()->SetRange(0, 0);
1541   meas->GetYaxis()->SetRange(0, 0);
1542
1543   gene->GetXaxis()->SetRangeUser(-10, 10);
1544   meas->GetXaxis()->SetRangeUser(-10, 10);
1545   AliPWG0Helper::CreateDividedProjections(gene, meas, "zy");
1546   gene->GetXaxis()->SetRange(0, 0);
1547   meas->GetXaxis()->SetRange(0, 0);
1548 }
1549
1550 TCanvas* Track2Particle2D(const char* fileName = "correction_map.root", const char* folder = "dndeta_correction")
1551 {
1552   gSystem->Load("libPWG0base");
1553
1554   Track2Particle2DCreatePlots(fileName);
1555
1556   TH2* corrYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx"));
1557   TH2* corrZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx"));
1558   TH2* corrZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy"));
1559
1560   Prepare2DPlot(corrYX);
1561   Prepare2DPlot(corrZX);
1562   Prepare2DPlot(corrZY);
1563
1564   const char* title = "";
1565   corrYX->SetTitle(title);
1566   corrZX->SetTitle(title);
1567   corrZY->SetTitle(title);
1568
1569   TCanvas* canvas = new TCanvas("Track2Particle2D", "Track2Particle2D", 1200, 400);
1570   canvas->Divide(3, 1);
1571
1572   canvas->cd(1);
1573   InitPadCOLZ();
1574   corrYX->Draw("COLZ");
1575
1576   canvas->cd(2);
1577   InitPadCOLZ();
1578   corrZX->Draw("COLZ");
1579
1580   canvas->cd(3);
1581   InitPadCOLZ();
1582   corrZY->Draw("COLZ");
1583
1584   canvas->SaveAs(Form("corr_track2particle_%d.gif", gMax));
1585   canvas->SaveAs(Form("corr_track2particle_%d.eps", gMax));
1586
1587   return canvas;
1588 }
1589
1590 void CompareTrack2Particle2D()
1591 {
1592   gSystem->Load("libPWG0base");
1593
1594   Track2Particle2DCreatePlots("correction_maponly-positive.root");
1595
1596   TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("pos_yx"));
1597   TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("pos_zx"));
1598   TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("pos_zy"));
1599
1600   Track2Particle2DCreatePlots("correction_maponly-negative.root");
1601
1602   TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("neg_yx"));
1603   TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("neg_zx"));
1604   TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("neg_zy"));
1605
1606   posYX->Divide(negYX);
1607   posZX->Divide(negZX);
1608   posZY->Divide(negZY);
1609
1610   Prepare2DPlot(posYX);
1611   Prepare2DPlot(posZX);
1612   Prepare2DPlot(posZY);
1613
1614   Float_t min = 0.8;
1615   Float_t max = 1.2;
1616
1617   posYX->SetMinimum(min);
1618   posYX->SetMaximum(max);
1619   posZX->SetMinimum(min);
1620   posZX->SetMaximum(max);
1621   posZY->SetMinimum(min);
1622   posZY->SetMaximum(max);
1623
1624   TCanvas* canvas = new TCanvas("CompareTrack2Particle2D", "CompareTrack2Particle2D", 1200, 400);
1625   canvas->Divide(3, 1);
1626
1627   canvas->cd(1);
1628   InitPadCOLZ();
1629   posYX->Draw("COLZ");
1630
1631   canvas->cd(2);
1632   InitPadCOLZ();
1633   posZX->Draw("COLZ");
1634
1635   canvas->cd(3);
1636   InitPadCOLZ();
1637   posZY->Draw("COLZ");
1638
1639   canvas->SaveAs("CompareTrack2Particle2D.gif");
1640   canvas->SaveAs("CompareTrack2Particle2D.eps");
1641 }
1642
1643 void Track2Particle3D()
1644 {
1645   // get left margin proper
1646
1647   TFile* file = TFile::Open("correction_map.root");
1648
1649   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1650
1651   corr->SetTitle("Correction Factor");
1652   SetRanges(corr->GetZaxis());
1653
1654   Prepare3DPlot(corr);
1655
1656   TCanvas* canvas = new TCanvas("Track2Particle3D", "Track2Particle3D", 500, 500);
1657   canvas->SetTheta(29.428);
1658   canvas->SetPhi(16.5726);
1659
1660   corr->Draw();
1661
1662   canvas->SaveAs("Track2Particle3D.gif");
1663   canvas->SaveAs("Track2Particle3D.eps");
1664 }
1665
1666 void Track2Particle3DAll()
1667 {
1668   TFile* file = TFile::Open("correction_map.root");
1669
1670   TH3* gene = dynamic_cast<TH3*> (file->Get("dndeta_correction/gene_nTrackToNPart"));
1671   TH3* meas = dynamic_cast<TH3*> (file->Get("dndeta_correction/meas_nTrackToNPart"));
1672   TH3* corr = dynamic_cast<TH3*> (file->Get("dndeta_correction/corr_nTrackToNPart"));
1673
1674   gene->SetTitle("Generated Particles");
1675   meas->SetTitle("Measured Tracks");
1676   corr->SetTitle("Correction Factor");
1677
1678   Prepare3DPlot(gene);
1679   Prepare3DPlot(meas);
1680   Prepare3DPlot(corr);
1681
1682   TCanvas* canvas = new TCanvas("Track2Particle3DAll", "Track2Particle3DAll", 1200, 400);
1683   canvas->Divide(3, 1);
1684
1685   canvas->cd(1);
1686   InitPad();
1687   gene->Draw();
1688
1689   canvas->cd(2);
1690   meas->Draw();
1691
1692   canvas->cd(3);
1693   corr->Draw();
1694
1695   canvas->SaveAs("Track2Particle3DAll.gif");
1696   canvas->SaveAs("Track2Particle3DAll.eps");
1697 }
1698
1699 void MultiplicityMC(Int_t xRangeMax = 50)
1700 {
1701   TFile* file = TFile::Open("multiplicityMC.root");
1702
1703   if (!file)
1704   {
1705     printf("multiplicityMC.root could not be opened.\n");
1706     return;
1707   }
1708
1709   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicityESD"));
1710   TH1F* fMultiplicityMC = dynamic_cast<TH1F*> (file->Get("fMultiplicityMC"));
1711   TH2F* fCorrelation = dynamic_cast<TH2F*> (file->Get("fCorrelation"));
1712
1713   TH1F* correction = new TH1F("MultiplicityMC_correction", "MultiplicityMC_correction;Ntracks;Npart", 76, -0.5, 75.5);
1714   TH1F* correctionWidth = new TH1F("MultiplicityMC_correctionwidth", "MultiplicityMC_correctionwidth;Ntracks;Npart", 76, -0.5, 75.5);
1715   //fMultiplicityMC->GetNbinsX(), fMultiplicityMC->GetXaxis()->GetXmin(), fMultiplicityMC->GetXaxis()->GetXmax());
1716   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1717   {
1718     TH1D* proj = fCorrelation->ProjectionX("_px", i, i+1);
1719     proj->Fit("gaus", "0");
1720     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1721     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1722
1723     continue;
1724
1725     // draw for debugging
1726     new TCanvas;
1727     proj->DrawCopy();
1728     proj->GetFunction("gaus")->DrawCopy("SAME");
1729   }
1730
1731   TH1F* fMultiplicityESDCorrected = new TH1F("fMultiplicityESDCorrected", "fMultiplicityESDCorrected", 2010, -0.5, 200.5);
1732
1733   for (Int_t i=1; i<=correction->GetNbinsX(); ++i)
1734   {
1735     Float_t mean = correction->GetBinContent(i);
1736     Float_t width = correctionWidth->GetBinContent(i);
1737
1738     Int_t fillBegin = fMultiplicityESDCorrected->FindBin(mean - width * 3);
1739     Int_t fillEnd   = fMultiplicityESDCorrected->FindBin(mean + width * 3);
1740     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1741
1742     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1743     {
1744       fMultiplicityESDCorrected->AddBinContent(j, TMath::Gaus(fMultiplicityESDCorrected->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fMultiplicityESD->GetBinContent(i));
1745     }
1746   }
1747
1748   TH1F* fMultiplicityESDCorrectedRebinned = dynamic_cast<TH1F*> (fMultiplicityESDCorrected->Clone("fMultiplicityESDCorrectedRebinned"));
1749   fMultiplicityESDCorrectedRebinned->Rebin(10);
1750   fMultiplicityESDCorrectedRebinned->Scale(0.1);
1751
1752   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
1753   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
1754   ratio->Divide(fMultiplicityMC);
1755
1756   TH1F* ratio2 = dynamic_cast<TH1F*> (fMultiplicityESDCorrectedRebinned->Clone("multiplicity_ratio_corrected"));
1757   ratio2->Divide(fMultiplicityMC);
1758
1759   TCanvas* canvas = new TCanvas("MultiplicityMC", "MultiplicityMC", 1500, 1000);
1760   canvas->Divide(3, 2);
1761
1762   fMultiplicityESD->GetXaxis()->SetRangeUser(0, xRangeMax);
1763   ratio->GetXaxis()->SetRangeUser(0, xRangeMax);
1764   fCorrelation->GetXaxis()->SetRangeUser(0, xRangeMax);
1765   fCorrelation->GetYaxis()->SetRangeUser(0, xRangeMax);
1766   correction->GetXaxis()->SetRangeUser(0, xRangeMax);
1767   fMultiplicityESDCorrected->GetXaxis()->SetRangeUser(0, xRangeMax);
1768   fMultiplicityESDCorrectedRebinned->GetXaxis()->SetRangeUser(0, xRangeMax);
1769
1770   canvas->cd(1); //InitPad();
1771   fMultiplicityESD->Draw();
1772   fMultiplicityMC->SetLineColor(2);
1773   fMultiplicityMC->Draw("SAME");
1774
1775   TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1776   legend->AddEntry(fMultiplicityESD, "ESD");
1777   legend->AddEntry(fMultiplicityMC, "MC");
1778   legend->Draw();
1779
1780   canvas->cd(2);
1781   fCorrelation->Draw("COLZ");
1782
1783   canvas->cd(3);
1784   correction->Draw();
1785   //correction->Fit("pol1");
1786   correctionWidth->SetLineColor(2);
1787   correctionWidth->Draw("SAME");
1788
1789   legend = new TLegend(0.2, 0.7, 0.45, 0.85);
1790   legend->AddEntry(correction, "#bar{x}");
1791   legend->AddEntry(correctionWidth, "#sigma");
1792   legend->Draw();
1793
1794   canvas->cd(4);
1795   ratio->Draw();
1796
1797   ratio2->SetLineColor(2);
1798   ratio2->Draw("SAME");
1799
1800   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1801   legend->AddEntry(ratio, "uncorrected");
1802   legend->AddEntry(ratio2, "corrected");
1803   legend->Draw();
1804
1805   canvas->cd(5);
1806   fMultiplicityESDCorrected->SetLineColor(kBlue);
1807   fMultiplicityESDCorrected->Draw();
1808   fMultiplicityMC->Draw("SAME");
1809   fMultiplicityESD->Draw("SAME");
1810
1811   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1812   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1813   legend->AddEntry(fMultiplicityMC, "MC");
1814   legend->AddEntry(fMultiplicityESD, "ESD");
1815   legend->Draw();
1816
1817   canvas->cd(6);
1818   fMultiplicityESDCorrectedRebinned->SetLineColor(kBlue);
1819   fMultiplicityESDCorrectedRebinned->Draw();
1820   fMultiplicityMC->Draw("SAME");
1821
1822   legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1823   legend->AddEntry(fMultiplicityESDCorrectedRebinned, "ESD corrected");
1824   legend->AddEntry(fMultiplicityMC, "MC");
1825   legend->Draw();
1826
1827   canvas->SaveAs("MultiplicityMC.gif");
1828 }
1829
1830 void MultiplicityESD()
1831 {
1832   TFile* file = TFile::Open("multiplicityESD.root");
1833
1834   if (!file)
1835   {
1836     printf("multiplicityESD.root could not be opened.\n");
1837     return;
1838   }
1839
1840   TH1F* fMultiplicityESD = dynamic_cast<TH1F*> (file->Get("fMultiplicity"));
1841
1842   TCanvas* canvas = new TCanvas("MultiplicityESD", "MultiplicityESD", 500, 500);
1843
1844   fMultiplicityESD->Draw();
1845 }
1846
1847 void CompareCorrection2Measured(Float_t ptMin = 0.301, const char* dataInput = "analysis_esd_raw.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1848 {
1849   loadlibs();
1850
1851   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1852   TFile::Open(correctionMapFile);
1853   dNdEtaCorrection->LoadHistograms();
1854
1855   TFile* file = TFile::Open(dataInput);
1856
1857   if (!file)
1858   {
1859     cout << "Error. File not found" << endl;
1860     return;
1861   }
1862
1863   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1864   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1865
1866   gROOT->cd();
1867   
1868   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
1869   hist1->SetTitle("mc");
1870   Printf("mc contains %f entries", hist1->Integral());
1871   Printf("mc contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), hist1->GetYaxis()->FindBin(-0.99), hist1->GetYaxis()->FindBin(0.99), hist1->GetZaxis()->FindBin(ptMin), hist1->GetNbinsZ()));
1872
1873   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
1874   hist2->SetTitle("esd");
1875   Printf("esd contains %f entries", hist2->Integral());
1876   Printf("esd contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), hist2->GetYaxis()->FindBin(-0.99), hist2->GetYaxis()->FindBin(0.99), hist2->GetZaxis()->FindBin(ptMin), hist2->GetNbinsZ()));
1877
1878   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1879   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
1880
1881   hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1882   hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1883   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
1884
1885   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
1886   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
1887   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
1888   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
1889   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
1890
1891   TH2* hist3 = (TH2*) dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram()->Clone("mc2");
1892   hist3->SetTitle("mc2");
1893   Printf("mc event contains %f entries", hist3->Integral());
1894   Printf("mc event contains %f entries in |vtx-z| < 10", hist3->Integral(hist3->GetXaxis()->FindBin(-9.9), hist3->GetXaxis()->FindBin(9.9), 1, hist3->GetNbinsY()));
1895
1896   TH2* hist4 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("esd2");
1897   hist4->SetTitle("esd2");
1898   Printf("esd event contains %f entries", hist4->Integral());
1899   Printf("esd event contains %f entries in |vtx-z| < 10", hist4->Integral(hist4->GetXaxis()->FindBin(-9.9), hist4->GetXaxis()->FindBin(9.9), 1, hist4->GetNbinsY()));
1900   
1901   ratio = (TH2*) hist3->Clone("ratio");
1902   ratio->Divide(hist4);
1903   
1904   new TCanvas; ratio->Draw("COLZ");
1905 }
1906
1907 void CompareCorrection2Generated(Float_t ptMin = 0.301, const char* dataInput = "analysis_mc.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1908 {
1909   loadlibs();
1910
1911   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1912   TFile::Open(correctionMapFile);
1913   dNdEtaCorrection->LoadHistograms();
1914
1915   TFile* file = TFile::Open(dataInput);
1916
1917   if (!file)
1918   {
1919     cout << "Error. File not found" << endl;
1920     return;
1921   }
1922
1923   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
1924   fdNdEtaAnalysis->LoadHistograms("dndetaTrVtx");
1925
1926   gROOT->cd();
1927   
1928   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram()->Clone("mc");
1929   hist1->SetTitle("mc");
1930   Printf("mc contains %f entries", hist1->Integral());
1931   Printf("mc contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), hist1->GetYaxis()->FindBin(-0.99), hist1->GetYaxis()->FindBin(0.99), hist1->GetZaxis()->FindBin(ptMin), hist1->GetNbinsZ()));
1932
1933   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram()->Clone("esd");
1934   hist2->SetTitle("esd");
1935   Printf("esd contains %f entries", hist2->Integral());
1936   Printf("esd contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), hist2->GetYaxis()->FindBin(-0.99), hist2->GetYaxis()->FindBin(0.99), hist2->GetZaxis()->FindBin(ptMin), hist2->GetNbinsZ()));
1937
1938   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1939   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "x");
1940
1941   hist1->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1942   hist2->GetXaxis()->SetRange(hist1->GetXaxis()->FindBin(-10), hist2->GetXaxis()->FindBin(10));
1943   AliPWG0Helper::CreateDividedProjections(hist1, hist2, "y");
1944
1945   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
1946   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
1947   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
1948   new TCanvas; gROOT->FindObject("mc_x_div_esd_x")->Draw("COLZ");
1949   new TCanvas; gROOT->FindObject("mc_y_div_esd_y")->Draw("COLZ");
1950 }
1951
1952 void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
1953 {
1954   loadlibs();
1955
1956   TFile* file = TFile::Open(dataInput);
1957
1958   if (!file)
1959   {
1960     cout << "Error. File not found" << endl;
1961     return;
1962   }
1963
1964   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
1965   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
1966
1967   TFile* file = TFile::Open(dataInput2);
1968
1969   if (!file)
1970   {
1971     cout << "Error. File not found" << endl;
1972     return;
1973   }
1974
1975   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
1976   fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
1977
1978   gROOT->cd();
1979
1980   TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
1981   hist1->SetTitle("esd1");
1982   Printf("esd1 contains %f entries", hist1->GetEntries());
1983   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()));
1984
1985   TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
1986   hist2->SetTitle("esd2");
1987   Printf("esd2 contains %f entries", hist2->GetEntries());
1988   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()));
1989
1990   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
1991
1992   new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
1993   new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
1994   new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
1995
1996   TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
1997   TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
1998
1999   Printf("event1 contains %f entries", event1->GetEntries());
2000   Printf("event2 contains %f entries", event2->GetEntries());
2001   Printf("event1 integral is %f", event1->Integral());
2002   Printf("event2 integral is %f", event2->Integral());
2003   Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
2004   Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
2005
2006   projx1 = event1->ProjectionX();
2007   projx2 = event2->ProjectionX();
2008
2009   new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
2010
2011   projx1->Divide(projx2);
2012   new TCanvas; projx1->Draw();
2013
2014   event1->Divide(event2);
2015   new TCanvas; event1->Draw("COLZ");
2016
2017 }
2018
2019 void DrawTrackletOrigin(const char* fileName = "correction_map.root", Bool_t myFile = kTRUE)
2020 {
2021   TFile::Open(fileName);
2022
2023   Int_t maxHists = 8;
2024   TH1* hist[8];
2025   
2026   const Int_t kRebin = 8;
2027
2028   const char* titles[] = { "PP", "SS", "PP'", "PS'", "PS", "SP'", "SS'", "" };
2029
2030   if (myFile)
2031   {
2032     for (Int_t i=0; i<maxHists; i++)
2033     {
2034       hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
2035       if (hist[i]->GetDimension() == 2)
2036         hist[i] = ((TH2*) hist[i])->ProjectionX(Form("fDeltaPhi_clone_%d", i));
2037     }
2038   }
2039   else
2040   {
2041     maxHists = 6;
2042     const char* names[] = { "DePhiPPTracklets", "DePhiSecTracklets", "DePhiPpTracklets", "DePhiPSTracklets", "DePhiPSdaugTracklets", "DePhiSPTracklets" }; 
2043     for (Int_t i=0; i<maxHists; i++)
2044       hist[i] = (TH1*) gFile->Get(names[i]);
2045   }
2046   
2047   // clone before rebinning
2048   good = (TH1*) hist[0]->Clone("good");
2049   good->Add(hist[4]);
2050   
2051   bad = (TH1*) hist[1]->Clone("bad");
2052   bad->Add(hist[2]);
2053   bad->Add(hist[3]);
2054   bad->Add(hist[5]);
2055   if (myFile)
2056     bad->Add(hist[6]);
2057   
2058   c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c");
2059   TH1* ref = 0;
2060   Bool_t nw = kFALSE;
2061   if (!c)
2062   {
2063     c = new TCanvas("c", "c", 600, 600);
2064     nw = kTRUE;
2065     ref = (TH1*) c->GetListOfPrimitives()->At(1);
2066   }  
2067   c->cd();
2068   c->SetRightMargin(0.05);
2069   c->SetTopMargin(0.05);
2070   c->SetLogy();
2071   c->SetGridx();
2072   c->SetGridy();
2073   
2074   Int_t order[] = { 0, 4, 1, 2, 3, 5, 6, 7 };
2075   //Int_t colors[]  = {1,2,4,1,2,4,1,2,4};
2076   Int_t colors[]  = {1,2,3,4,6,7,8,102};
2077   Int_t markers[]  = {20, 21, 22, 23, 24, 25, 26, 27, 28};
2078   
2079   TLegend* legend = new TLegend(0.75, 0.6, 0.93, 0.93);
2080   legend->SetFillColor(0);
2081   legend->SetTextSize(0.04);
2082
2083   Int_t total = 0;
2084   for (Int_t ii=0; ii<maxHists; ii++)
2085   {
2086     i = order[ii];
2087     
2088     hist[i]->Rebin(kRebin);
2089     hist[i]->SetStats(kFALSE);
2090     hist[i]->SetLineColor(colors[i]);
2091     hist[i]->SetLineWidth(2);
2092     //hist[i]->SetMarkerStyle(markers[i]);
2093     //hist[i]->SetMarkerColor(colors[i]);
2094     //hist[i]->SetLineStyle(ii+1);
2095     hist[i]->GetXaxis()->SetRangeUser(-0.09, 0.09);
2096     hist[i]->GetYaxis()->SetRangeUser(5, hist[i]->GetMaximum() * 2);
2097     hist[i]->GetYaxis()->SetTitleOffset(1.3);
2098     hist[i]->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2099     
2100     if (i == 0 && ref)
2101       hist[i]->Scale(1.0 / hist[i]->GetMaximum() * ref->GetMaximum());
2102     
2103     hist[i]->DrawCopy(((i == 0 && nw) ? "" : "SAME"));
2104
2105     total += hist[i]->GetEntries();
2106
2107     if (i != 7)
2108       legend->AddEntry(hist[i], titles[i], "L");
2109   }
2110
2111   legend->Draw();
2112   c->SaveAs("spd_tracklets_deltaphi_detailed.eps");
2113
2114   Printf("Total: %d", total);
2115   for (Int_t i=0; i<maxHists; i++)
2116     Printf("Histogram %d (%s) contains %.2f %% of the entries", i, titles[i], 100.0 * hist[i]->GetEntries() / total);
2117
2118   printf("|  Delta phi  |  Acc. %%  |  ");
2119   for (Int_t i=0; i<maxHists; i++)
2120     printf("%3s %%   |  ", titles[i]);
2121   Printf("");
2122
2123   for (Float_t f = 0.01; f < 0.09; f += 0.01)
2124   {
2125     Int_t integralBegin = hist[0]->GetXaxis()->FindBin(-f);
2126     Int_t integralEnd = hist[0]->GetXaxis()->FindBin(f);
2127
2128     Int_t total2 = 0;
2129     for (Int_t i=0; i<maxHists; i++)
2130       total2 += (Int_t) hist[i]->Integral(integralBegin, integralEnd);
2131
2132     printf("|    %.2f     |  %6.2f  |  ", f, 100.0 * total2 / total);
2133
2134     for (Int_t i=0; i<maxHists; i++)
2135       printf("%6.2f  |  ", (hist[i]->GetEntries() > 0) ? (100.0 * hist[i]->Integral(integralBegin, integralEnd) / hist[i]->GetEntries()) : -1.0);
2136     Printf("");
2137   }
2138   
2139   eff = new TH1F("eff", ";#Delta#varphi cut (rad.)", 101,-0.0005, 0.1005);
2140   cont = new TH1F("cont", "cont", 101,-0.0005, 0.1005);
2141   signalOverBg = new TH1F("signalOverBg", "signalOverBg", 101,-0.0005, 0.1005);
2142   for (Float_t cut=0.000; cut<0.10; cut += 0.001)
2143   {
2144     Float_t accGood = good->Integral(good->GetXaxis()->FindBin(-cut), good->GetXaxis()->FindBin(cut));
2145     Float_t accBad = bad->Integral(bad->GetXaxis()->FindBin(-cut), bad->GetXaxis()->FindBin(cut));
2146     Float_t sB = accGood / accBad;
2147     eff->Fill(cut, 100.0 * accGood / good->Integral());
2148     cont->Fill(cut, 100.0 * accBad / (accGood + accBad));
2149     signalOverBg->Fill(cut, sB);
2150   }
2151   
2152   //new TCanvas; signalOverBg->Draw();
2153   
2154   c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c2");
2155   Bool_t nw = kFALSE;
2156   if (!c)
2157   {
2158     c = new TCanvas("c2", "c2", 600, 600);
2159     nw = kTRUE;
2160   }
2161   c->cd();
2162   c->SetRightMargin(0.05);
2163   c->SetTopMargin(0.05);
2164   c->SetGridx();
2165   c->SetGridy();
2166   gPad->SetLogy();
2167   good->Rebin(kRebin);
2168   bad->Rebin(kRebin);
2169   good->GetXaxis()->SetRangeUser(-0.09, 0.09);
2170   good->GetYaxis()->SetTitleOffset(1.3);
2171   good->SetStats(0);
2172   good->GetXaxis()->SetTitle("#Delta#varphi (rad.)");  
2173   good->DrawCopy((nw) ? "" : "SAME");
2174   
2175   bad->SetLineColor(2);
2176   bad->SetLineStyle(2);
2177   bad->SetLineWidth(2);
2178   //bad->SetMarkerColor(2);
2179   //bad->SetMarkerStyle(7);
2180   bad->DrawCopy("SAME");
2181   
2182   TLegend* legend = new TLegend(0.2, 0.13, 0.85, 0.25);
2183   legend->SetFillColor(0);
2184   legend->SetTextSize(0.04);
2185   legend->AddEntry(good, "Primaries", "L");
2186   legend->AddEntry(bad, "Secondaries + Background", "L");
2187   legend->Draw();
2188   
2189   c->SaveAs("spd_tracklets_deltaphi.eps");
2190   
2191   c = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c3");
2192   Bool_t nw = kFALSE;
2193   if (!c)
2194   {
2195     c = new TCanvas("c3", "c3", 600, 600);
2196     nw = kTRUE;
2197   }
2198   c->cd();
2199   c->SetRightMargin(0.05);
2200   c->SetTopMargin(0.05);
2201   c->SetGridx();
2202   c->SetGridy();
2203   
2204   TLegend* legend = new TLegend(0.5, 0.6, 0.93, 0.75);
2205   legend->SetFillColor(0);
2206   legend->SetTextSize(0.04);
2207   legend->AddEntry(eff, "Efficiency (%)", "L");
2208   legend->AddEntry(cont, "Contamination (%)", "L");
2209   
2210   eff->SetStats(0);
2211   eff->GetXaxis()->SetRangeUser(0, 0.08);
2212   eff->GetYaxis()->SetRangeUser(1e-3, 105);
2213   eff->SetLineWidth(2);
2214   eff->DrawCopy((nw) ? "" : "SAME");
2215   cont->SetLineStyle(2);
2216   cont->SetLineWidth(2);
2217   cont->SetLineColor(2);
2218   cont->DrawCopy("SAME");
2219   legend->Draw();
2220   
2221   c->SaveAs("spd_tracklets_efficiency.eps");
2222 }
2223
2224 void DrawTrackletOrigin_Compare(const char* file1, const char* file2)
2225 {
2226   DrawTrackletOrigin(file1);
2227   good1 = (TH1*) gROOT->FindObject("good")->Clone("good1");
2228   bad1 = (TH1*) gROOT->FindObject("bad")->Clone("bad1");
2229
2230   DrawTrackletOrigin(file2);
2231   good2 = (TH1*) gROOT->FindObject("good")->Clone("good2");
2232   bad2 = (TH1*) gROOT->FindObject("bad")->Clone("bad2");
2233      
2234   c = new TCanvas("c4", "c4", 600, 600);
2235   c->SetRightMargin(0.05);
2236   c->SetTopMargin(0.05);
2237   c->SetGridx();
2238   c->SetGridy();
2239   gPad->SetLogy();
2240   
2241   good1->Draw();
2242   bad1->SetLineColor(1);
2243   bad1->SetMarkerColor(1);
2244   bad1->Draw("SAME");
2245   
2246   Float_t factor = (good1->Integral() + bad1->Integral()) / (good2->Integral() + bad2->Integral());
2247   
2248   good2->Scale(factor);
2249   bad2->Scale(factor);
2250   
2251   good2->SetLineColor(2);
2252   bad2->SetMarkerColor(2);
2253   
2254   good2->Draw("SAME");
2255   bad2->Draw("SAME");
2256   
2257   good1->GetYaxis()->SetRangeUser(1, TMath::Max(good1->GetMaximum(), good2->GetMaximum()) * 1.1);
2258 }
2259   
2260 void Tracklets_Asymmetry()
2261 {
2262   TFile::Open("correction_map.root");
2263
2264   Int_t maxHists = 7;
2265   TH1* hist[8];
2266
2267   Int_t colors[]  = {1,2,3,4,6,7,8,102};
2268   const char* titles[] = { "PP", "SS", "PP'", "PS'", "PS", "SP'", "SS'", "" };
2269
2270   TLegend* legend = new TLegend(0.75, 0.6, 0.93, 0.93);
2271   
2272   for (Int_t i=0; i<maxHists; i++)
2273   {
2274     hist[i] = (TH1*) gFile->Get(Form("fDeltaPhi_%d", i));
2275     hist[i]->Rebin(10);
2276     
2277     for (Int_t j=hist[i]->GetNbinsX()/2; j<=hist[i]->GetNbinsX(); j++)
2278       if (hist[i]->GetBinContent(j) > 0)
2279         hist[i]->SetBinContent(j, (hist[i]->GetBinContent(j) -  hist[i]->GetBinContent(hist[i]->GetXaxis()->FindBin(-hist[i]->GetXaxis()->GetBinCenter(j)))) / hist[i]->GetBinContent(j));
2280       
2281     hist[i]->SetStats(kFALSE);
2282     hist[i]->SetLineColor(colors[i]);
2283     hist[i]->GetXaxis()->SetRangeUser(0.001, 0.09);
2284     //hist[i]->GetYaxis()->SetRangeUser(5, hist[i]->GetMaximum() * 2);
2285     hist[i]->GetYaxis()->SetTitleOffset(1.3);
2286     hist[i]->GetXaxis()->SetTitle("#Delta#varphi (rad.)");
2287     hist[i]->Draw(((i == 0) ? "" : "SAME"));
2288     
2289     legend->AddEntry(hist[i], titles[i], "L");
2290   }
2291   
2292   legend->Draw();
2293 }
2294
2295 TH2* GetCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2296 {
2297   // returns the correction factor with pt integrated out
2298
2299   loadlibs();
2300
2301   TFile::Open(fileName);
2302
2303   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
2304   if (!dNdEtaCorrection->LoadHistograms())
2305     return;
2306
2307   //  hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
2308
2309   gener = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
2310   measu = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
2311
2312   gener->GetZaxis()->SetRange(gener->GetZaxis()->FindBin(ptmin), gener->GetNbinsZ()+1);
2313   TH2D *gener_xy = gener->Project3D("yx");
2314
2315   measu->GetZaxis()->SetRange(measu->GetZaxis()->FindBin(ptmin), measu->GetNbinsZ()+1);
2316   TH2D *measu_xy = measu->Project3D("yx");
2317
2318   cout << measu->GetZaxis()->FindBin(ptmin) << " " << measu->GetNbinsZ()+1 << endl;
2319
2320   TCanvas *canp = new TCanvas("canp","canp",600,1000);
2321   canp->Divide(1,2,0.0001,0.0001);
2322   canp->cd(1);
2323   gener_xy->Draw("COLZ");
2324   canp->cd(2);
2325   measu_xy->Draw("COLZ");
2326
2327
2328   TCanvas *canpr = new TCanvas("canpr","canpr",700,500);
2329   canpr->cd();
2330   TH2D *proj = new TH2D(*gener_xy);
2331   proj->Divide(measu_xy);
2332
2333 //   proj = hist->Project3D("yx");
2334   proj->Draw("COLZ");
2335
2336   return proj;
2337 }
2338
2339 void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2340 {
2341   TH2* proj = GetCorrection(fileName, dirName, ptmin);
2342
2343   const Float_t limit = 5;
2344
2345   TString array = "{";
2346   TString arrayEnd = "}";
2347
2348   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2349   {
2350     Int_t begin = -1;
2351     Int_t end = -1;
2352     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2353     {
2354       if (begin == -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2355         begin = x;
2356       if (begin != -1 && proj->GetBinContent(x, y) > 0 && proj->GetBinContent(x, y) < limit)
2357         end = x;
2358     }
2359     Printf("Limits for y = %d are %d to %d", y, begin, end);
2360
2361     if (y > 1)
2362       array += ", ";
2363     array += Form("%d", begin);
2364
2365     if (y > 1)
2366       arrayEnd.Prepend(", ");
2367     arrayEnd.Prepend(Form("%d", (end == -1) ? -1 : proj->GetNbinsX() + 1 - end));
2368   }
2369   array += "}";
2370   arrayEnd.Prepend("{");
2371
2372   Printf("Begin array:");
2373   Printf("%s", array.Data());
2374
2375   Printf("End array (mirrored) (should be the same):");
2376   Printf("%s", arrayEnd.Data());
2377 }
2378
2379 void AverageMultiplicity(const char* fileName = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
2380 {
2381   loadlibs();
2382
2383   TFile::Open(fileName);
2384
2385   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
2386   dNdEtaCorrection->LoadHistograms();
2387   TH2* events = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram();
2388   TH3* tracks = dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram();
2389
2390   Float_t nEvents = events->Integral(events->GetXaxis()->FindBin(-1), events->GetXaxis()->FindBin(1), 0, events->GetNbinsY()+1);
2391   Float_t nTracks = tracks->Integral(tracks->GetXaxis()->FindBin(-1), tracks->GetXaxis()->FindBin(1), tracks->GetYaxis()->FindBin(-0.39), tracks->GetYaxis()->FindBin(0.59), 0, tracks->GetNbinsZ()+1);
2392
2393   Printf("%f %f --> %f", nEvents, nTracks, nTracks / nEvents);
2394 }
2395
2396 void GetAverageCorrectionFactor(Float_t etaRange = 1.5, Float_t vertexRange = 9.9, const char* rawFile = "analysis_esd_raw.root", const char* mcFile = "analysis_mc.root")
2397 {
2398   loadlibs();
2399
2400   TFile::Open(rawFile);
2401   dNdEtaAnalysis* raw = new dNdEtaAnalysis("dndeta", "dndeta");
2402   raw->LoadHistograms("fdNdEtaAnalysisESD");
2403   raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->GetXaxis()->SetRangeUser(-vertexRange, vertexRange);
2404   tracks = raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Project3D("y");
2405   events = raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("events", 0, raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->GetNbinsY() + 1);
2406   Float_t nEvents = events->Integral(events->FindBin(-vertexRange), events->FindBin(vertexRange));
2407   tracks->Scale(1.0 / nEvents / tracks->GetBinWidth(1));
2408
2409   TFile::Open(mcFile);
2410   dNdEtaAnalysis* mc = new dNdEtaAnalysis("dndeta", "dndeta");
2411   mc->LoadHistograms("dndetaTrVtx");
2412   mcH = mc->GetdNdEtaPtCutOffCorrectedHistogram(0);
2413
2414   new TCanvas;
2415   mcH->SetLineColor(2);
2416   mcH->DrawCopy();
2417   tracks->DrawCopy("SAME");
2418
2419   new TCanvas;
2420   mcH->GetYaxis()->SetRangeUser(0, 5);
2421   mcH->Divide(tracks);
2422   mcH->DrawCopy();
2423   mcH->Fit("pol0", "", "", -etaRange, etaRange);
2424 }
2425
2426 void TrackCuts_Comparison_MC(char* histName, Int_t plotWhich = 0, const char* fileName = "correction_map.root", Bool_t mirror = kFALSE)
2427 {
2428   // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
2429   //    --> manually disable it in the run.C
2430   //
2431   // plotWhich: 0 = only before
2432   //            1 = both
2433   //            2 = only after
2434   //
2435   // mirror: kTRUE --> project negative values on the positive side
2436   
2437
2438   file = TFile::Open(fileName);
2439
2440   Int_t count = 0;
2441   Int_t colors[] = { 1, 2, 3, 4, 5, 6 };
2442
2443   TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
2444   TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
2445   TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
2446
2447   TCanvas* c1 = new TCanvas(histName, histName, 800, 1200);
2448   c1->Divide(1, 2);
2449   //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
2450   //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
2451
2452   const char* folders2[] = { "before_cuts", "after_cuts" };
2453   Bool_t first = kTRUE;
2454   for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
2455   {
2456     const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
2457     const char* names[] =    { "all", "primaries", "secondaries" };
2458     TH1* base = 0;
2459     TH1* prim = 0;
2460     TH1* sec = 0;
2461     for (Int_t i = 0; i < 3; i++)
2462     {
2463       TString folder;
2464       folder.Form("%s/%s/%s", folders1[i], folders2[j], histName);
2465       TH1* hist = (TH1*) file->Get(folder);
2466       
2467       if (mirror)
2468       {
2469         for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
2470         {
2471           Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
2472           if (bin != newBin)
2473           {
2474             hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
2475             hist->SetBinContent(bin, 0);
2476           }
2477         }
2478       }
2479       
2480       legend->AddEntry(hist, Form("%s %s", names[i], folders2[j]));
2481
2482       c1->cd(1);
2483       hist->SetLineColor(colors[count]);
2484       hist->DrawCopy((count == 0) ? "" : "SAME");
2485
2486       switch (i)
2487       {
2488         case 0: base = hist; break;
2489         case 1: prim = hist; break;
2490         case 2: sec = hist; break;
2491       }
2492
2493       count++;
2494     }
2495
2496     TH1* eff    = (TH1*) prim->Clone("eff"); eff->Reset();
2497     TH1* purity = (TH1*) prim->Clone("purity"); purity->Reset();
2498
2499     for (Int_t bin = 1; bin <= prim->GetNbinsX(); bin++)
2500     {
2501       eff->SetBinContent(bin, prim->Integral(1, bin) / prim->Integral(1, prim->GetNbinsX() + 1));
2502       if (prim->Integral(1, bin) + sec->Integral(1, bin) > 0)
2503         purity->SetBinContent(bin, sec->Integral(1, bin) / (prim->Integral(1, bin) + sec->Integral(1, bin)));
2504     }
2505
2506     eff->GetYaxis()->SetRangeUser(0, 1);
2507     eff->SetLineColor(colors[0+j*2]);
2508     eff->SetStats(kFALSE);
2509     purity->SetLineColor(colors[1+j*2]);
2510
2511     legend3->AddEntry(eff, Form("%s: efficiency", folders2[j]));
2512     legend3->AddEntry(purity, Form("%s: contamination", folders2[j]));
2513
2514     c1->cd(2);
2515     eff->DrawCopy((first) ? "" : "SAME");
2516     first = kFALSE;
2517     purity->DrawCopy("SAME");
2518   }
2519
2520   c1->cd(1)->SetLogy();
2521   c1->cd(1)->SetGridx();
2522   c1->cd(1)->SetGridy();
2523   legend->Draw();
2524
2525   //c2->cd();
2526  // c2->SetGridx();
2527  // c2->SetGridy();
2528   //legend2->Draw();
2529
2530   c1->cd(2)->SetGridx();
2531   c1->cd(2)->SetGridy();
2532   legend3->Draw();
2533
2534   //c1->SaveAs(Form("%s.png", histName));
2535 }
2536
2537 void TrackCuts_Comparison_Data(char* histName, Int_t plotWhich, const char* fileName1, const char* fileName2, Bool_t mirror = kFALSE, const char* label1 = "file1", const char* label2 = "file2")
2538 {
2539   // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
2540   //    --> manually disable it in the run.C
2541   //
2542   // plotWhich: 0 = only before
2543   //            1 = both
2544   //            2 = only after
2545   //
2546   // mirror: kTRUE --> project negative values on the positive side
2547   
2548
2549   Int_t count = 0;
2550   Int_t colors[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
2551
2552   TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
2553   legend->SetFillColor(0);
2554   TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
2555   TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
2556
2557   TCanvas* c1 = new TCanvas(histName, histName, 600, 600);
2558   //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
2559   //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
2560
2561   const char* folders2[] = { "before_cuts", "after_cuts" };
2562   Bool_t first = kTRUE;
2563   for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
2564   {
2565     const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
2566     const char* names[] =    { "all", "primaries", "secondaries" };
2567     
2568     Float_t normalize[3];
2569     
2570     for (Int_t i = 0; i < 2; i++)
2571     {
2572       file = TFile::Open((i == 0) ? fileName1 : fileName2);
2573       
2574       for (Int_t k = 1; k < 3; k++)
2575       {
2576         TString folder;
2577         folder.Form("%s/%s/%s", folders1[k], folders2[j], histName);
2578         Printf("%s", folder.Data());
2579         TH1* hist = (TH1*) file->Get(folder);
2580         
2581         if (mirror)
2582         {
2583           for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
2584           {
2585             Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
2586             if (bin != newBin)
2587             {
2588               hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
2589               hist->SetBinContent(bin, 0);
2590             }
2591           }
2592         }
2593       
2594         if (i == 0)
2595         {
2596           normalize[k] = hist->Integral();
2597         }
2598         else
2599           hist->Scale(normalize[k] / hist->Integral());
2600         
2601         legend->AddEntry(hist, Form("%s %s %s", (i == 0) ? label1 : label2, (k == 1) ? "primaries" : "secondaries", folders2[j]));
2602   
2603         c1->cd();
2604         hist->SetStats(0);
2605         hist->SetLineColor(colors[count]);
2606         hist->DrawCopy((count == 0) ? "" : "SAME");
2607   
2608         count++;
2609       }
2610     }
2611
2612   }
2613
2614   //c1->SetLogy();
2615   c1->SetGridx();
2616   c1->SetGridy();
2617   legend->Draw();
2618 }
2619
2620 void TrackCuts_DCA()
2621 {
2622   file = TFile::Open("correction_map.root");
2623   hist = (TH2*) file->Get("esd_track_cuts/before_cuts/dXYvsDZ");
2624
2625   TCanvas* c1 = new TCanvas("c1", "c1", 600, 600);
2626   c1->SetLogz();
2627   c1->SetRightMargin(0.12);
2628   c1->SetBottomMargin(0.12);
2629
2630   hist->SetStats(kFALSE);
2631   hist->Draw("COLZ");
2632
2633   ellipse = new TEllipse(0, 0, 4);
2634   ellipse->SetLineWidth(2);
2635   ellipse->SetLineStyle(2);
2636   ellipse->SetFillStyle(0);
2637   ellipse->Draw();
2638
2639   c1->SaveAs("trackcuts_dca_2d.eps");
2640 }
2641
2642 void FindNSigma(TH2* hist, Int_t nSigma = 1)
2643 {
2644   TH1* proj = hist->ProjectionY();
2645   proj->Reset();
2646
2647   for (Int_t bin=1; bin<=proj->GetNbinsX(); bin++)
2648   {
2649     if (hist->Integral(1, hist->GetNbinsX(), bin, bin) == 0)
2650       continue;
2651
2652     Int_t limit = -1;
2653     for (limit = 1; limit<=hist->GetNbinsX(); limit++)
2654   }
2655 }
2656
2657 void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
2658 {
2659   TH2* proj = GetCorrection(fileName, dirName, ptmin);
2660
2661   for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
2662     for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
2663       if (proj->GetBinContent(x, y) > 5 || proj->GetBinContent(x, y) == 0)
2664       {
2665         proj->SetBinContent(x, y, 0);
2666       }
2667       else
2668         proj->SetBinContent(x, y, 1);
2669
2670
2671   input->Multiply(proj);
2672 }
2673
2674 void MakeGaussianProfile(const char* histName = "fVertexCorrelation", Bool_t subtractMean = kFALSE)
2675 {
2676     TFile::Open("correction_map.root");
2677
2678     TH2* hist2d = (TH2*) gFile->Get(histName);
2679     hist2d->Sumw2();
2680
2681     TH1* result = hist2d->ProjectionX("result");
2682     result->GetYaxis()->SetTitle(hist2d->GetYaxis()->GetTitle());
2683     result->Reset();
2684
2685     for (Int_t x=1; x<hist2d->GetNbinsX(); ++x)
2686     {
2687         hist = hist2d->ProjectionY(Form("temp_%d", x), x, x);
2688         if (hist->GetEntries() == 0)
2689             continue;
2690         if (hist->Fit("gaus") == 0)
2691         {
2692             func = hist->GetFunction("gaus");
2693             mean = func->GetParameter(1);
2694             error = func->GetParError(1);
2695
2696             if (subtractMean)
2697                 mean = hist2d->GetXaxis()->GetBinCenter(x) - mean;
2698
2699             result->SetBinContent(x, mean);
2700             result->SetBinError(x, error);
2701
2702             if (x % 10 == 0)
2703             {
2704                 new TCanvas;
2705                 ((TH1*) hist->Clone())->DrawCopy();
2706             }
2707         }
2708         //break;
2709     }
2710
2711     new TCanvas;
2712     result->GetYaxis()->SetRangeUser(-0.2, 0.2);
2713     result->Draw();
2714 }
2715
2716 TH2* GetAcceptance(void* corr2d_void)
2717 {
2718         corr2d = (AliCorrectionMatrix2D*) corr2d_void;
2719         corr_xy = (TH2*) corr2d->GetCorrectionHistogram()->Clone("acceptance");
2720
2721         // fold in acceptance
2722         for (Int_t x=1; x<=corr_xy->GetNbinsX(); ++x)
2723                 for (Int_t y=1; y<=corr_xy->GetNbinsY(); ++y)
2724                 {
2725                         if (corr_xy->GetBinContent(x, y) > 1.5)
2726                                 corr_xy->SetBinContent(x, y, 0);
2727
2728                         if (corr_xy->GetBinContent(x, y) > 0)
2729                                 corr_xy->SetBinContent(x, y, 1);
2730
2731                         corr_xy->SetBinError(x, y, 0);
2732                 }
2733
2734         return corr_xy;
2735 }
2736
2737 void ZeroOutsideAcceptance(TH2* acc, TH3* hist)
2738 {
2739   for (Int_t x=0; x<=acc->GetNbinsX()+1; ++x)
2740     for (Int_t y=0; y<=acc->GetNbinsY()+1; ++y)
2741     {
2742       if (acc->GetBinContent(x, y) > 1.5 || acc->GetBinContent(x, y) == 0)
2743       {
2744         for (Int_t z=0; z<=hist->GetNbinsZ()+1; ++z)
2745         {
2746           hist->SetBinContent(x, y, z, 0);
2747           hist->SetBinError(x, y, z, 0);
2748         }
2749       }
2750     }
2751 }
2752
2753 void DrawPhi()
2754 {
2755   loadlibs();
2756
2757   TFile::Open("correction_map.root");
2758   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
2759   if (!dNdEtaCorrection->LoadHistograms())
2760     return 0;
2761
2762   TFile::Open("analysis_esd.root");
2763   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
2764   fdNdEtaAnalysis->LoadHistograms();
2765
2766   // acc. map!
2767   //acc = GetAcceptance(dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000));
2768   acc = dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000)->GetCorrectionHistogram();
2769   //new TCanvas; acc->Draw("COLZ");
2770
2771   histG = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
2772   ZeroOutsideAcceptance(acc, histG);
2773   //new TCanvas; histG->Project3D("yx")->Draw("COLZ");
2774   //histG->GetYaxis()->SetRangeUser(-0.9, 0.9);
2775
2776   histM = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
2777   ZeroOutsideAcceptance(acc, histM);
2778   //histM->GetYaxis()->SetRangeUser(-0.9, 0.9);
2779
2780   TFile::Open("analysis_mc.root");
2781   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndetaTrVtxMC", "dndetaTrVtxMC");
2782   fdNdEtaAnalysis2->LoadHistograms("dndetaTrVtx");
2783
2784   histMC = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
2785   ZeroOutsideAcceptance(acc, histMC);
2786   //new TCanvas; histMC->Project3D("yx2")->Draw("COLZ");
2787
2788   //histG->GetZaxis()->SetRangeUser(1,2); histMC->GetZaxis()->SetRangeUser(1,2);
2789   new TCanvas; a = histG->Project3D("yx3"); a->Add(histMC->Project3D("yx4"), -1); a->Draw("COLZ");
2790
2791   //histMC->GetYaxis()->SetRangeUser(-0.9, 0.9);
2792
2793   c = new TCanvas;
2794
2795   histG->GetXaxis()->SetRangeUser(-9.9, 9.9);
2796   histG->Project3D("z")->Draw();
2797
2798   histM->GetXaxis()->SetRangeUser(-9.9, 9.9);
2799   proj = histM->Project3D("z2");
2800   proj->SetLineColor(2);
2801   proj->Draw("SAME");
2802
2803   histMC->GetXaxis()->SetRangeUser(-9.9, 9.9);
2804   projMC = histMC->Project3D("z3");
2805   projMC->SetLineColor(3);
2806   projMC->Draw("SAME");
2807 }
2808
2809 void PrintEventStats(Int_t corrID = 3)
2810 {
2811   loadlibs();
2812
2813   /*
2814   TFile::Open("analysis_mc.root");
2815   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD");
2816   fdNdEtaAnalysis->LoadHistograms();
2817   trackHist = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
2818   eventHist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetGeneratedHistogram();
2819   */
2820
2821   TFile::Open("correction_map.root");
2822   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
2823   if (!dNdEtaCorrection->LoadHistograms())
2824     return;
2825   trackHist = dNdEtaCorrection->GetCorrection(corrID)->GetTrackCorrection()->GetGeneratedHistogram();
2826   eventHist = dNdEtaCorrection->GetCorrection(corrID)->GetEventCorrection()->GetGeneratedHistogram();
2827
2828   trackHist->GetXaxis()->SetRange(0, trackHist->GetNbinsX()+1);
2829   trackHist->GetZaxis()->SetRange(0, trackHist->GetNbinsZ()+1);
2830   eta = trackHist->Project3D("y");
2831
2832   events = eventHist->Integral(0, eventHist->GetNbinsX()+1, 0, eventHist->GetNbinsY()+1);
2833
2834   eta->Scale(1.0 / events);
2835
2836   Float_t avgN = eta->Integral(0, eta->GetNbinsX()+1);
2837   Printf("<N> = %f", avgN);
2838
2839   eta->Scale(1.0 / eta->GetXaxis()->GetBinWidth(1));
2840
2841   Printf("dndeta | eta = 0 is %f", (eta->GetBinContent(eta->FindBin(0.01)) + eta->GetBinContent(eta->FindBin(-0.01))) / 2);
2842   Printf("dndeta in |eta| < 0.5 is %f", eta->Integral(eta->FindBin(-0.49), eta->FindBin(0.49)) / (eta->FindBin(0.49) - eta->FindBin(-0.49) + 1));
2843   Printf("dndeta in |eta| < 1 is %f", eta->Integral(eta->FindBin(-0.99), eta->FindBin(0.99)) / (eta->FindBin(0.99) - eta->FindBin(-0.99) + 1));
2844   Printf("dndeta in |eta| < 1.5 is %f", eta->Integral(eta->FindBin(-1.49), eta->FindBin(1.49)) / (eta->FindBin(1.49) - eta->FindBin(-1.49) + 1));
2845
2846   stats = (TH2*) gFile->Get("fEventStats");
2847   proj = stats->ProjectionX();
2848   gROOT->ProcessLine(".L PrintHist.C");
2849   PrintHist2D(stats);
2850   PrintHist(proj);
2851   
2852   Printf("+++ TRIGGER EFFICIENCIES +++");
2853   
2854   Printf("INEL = %.1f", 100. * (proj->GetBinContent(1) - stats->GetBinContent(1, 1) - stats->GetBinContent(1, 3)) / proj->GetBinContent(1));
2855   Printf("NSD  = %.1f", 100. * (proj->GetBinContent(2) - stats->GetBinContent(2, 1) - stats->GetBinContent(2, 3)) / proj->GetBinContent(2));
2856   Printf("ND  = %.1f",  100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1) - stats->GetBinContent(3, 3)) / proj->GetBinContent(3));
2857   Printf("SD  = %.1f",  100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1) - stats->GetBinContent(4, 3)) / proj->GetBinContent(4));
2858   Printf("DD  = %.1f",  100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1) - stats->GetBinContent(5, 3)) / proj->GetBinContent(5));
2859   
2860   Printf("+++ TRIGGER + VERTEX EFFICIENCIES +++");
2861   
2862   Printf("INEL = %.1f", 100. * stats->GetBinContent(1, 4) / proj->GetBinContent(1));
2863   Printf("NSD  = %.1f", 100. * stats->GetBinContent(2, 4) / proj->GetBinContent(2));
2864   Printf("ND  = %.1f",  100. * stats->GetBinContent(3, 4) / proj->GetBinContent(3));
2865   Printf("SD  = %.1f",  100. * stats->GetBinContent(4, 4) / proj->GetBinContent(4));
2866   Printf("DD  = %.1f",  100. * stats->GetBinContent(5, 4) / proj->GetBinContent(5));
2867   
2868   for (Int_t i=7; i<=proj->GetNbinsX(); i++)
2869     if (proj->GetBinContent(i) > 0)
2870       Printf("bin %d (process type %d) = %.2f", i, (Int_t) proj->GetXaxis()->GetBinCenter(i), 100.0 * (proj->GetBinContent(i) - stats->GetBinContent(i, 1)) / proj->GetBinContent(i));
2871   
2872   //eta->Draw();
2873 }
2874
2875 void TestAsymmetry()
2876 {
2877   loadlibs();
2878
2879   TFile* file2 = TFile::Open("analysis_mc.root");
2880   
2881   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
2882   fdNdEtaAnalysis->LoadHistograms();
2883   fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
2884   
2885   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy();
2886   
2887   hist = (TH1*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
2888   hist2 = (TH1*) hist->Clone("hist2");
2889   for (Int_t x=1; x<=hist->GetNbinsX(); x++)
2890     for (Int_t y=1; y<=hist->GetNbinsY(); y++)
2891       for (Int_t z=1; z<=hist->GetNbinsZ(); z++)
2892       {
2893         Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
2894         hist->SetBinContent(x, y, z, hist2->GetBinContent(hist->GetNbinsX() / 2, TMath::Min(y, hist->GetNbinsY() + 1 - y), z));
2895       }
2896   
2897   hist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
2898   for (Int_t x=1; x<=hist->GetNbinsX(); x++)
2899     for (Int_t y=1; y<=hist->GetNbinsY(); y++)
2900       {
2901         //Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
2902         hist->SetBinContent(x, y, hist->GetBinContent(hist->GetNbinsX() / 2, y));
2903       }
2904   
2905   fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
2906   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerColor(2);
2907   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetLineColor(2);
2908   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerStyle(5);
2909   fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy("SAMEP");
2910 }
2911
2912 void DeltaPhiFromPt(Float_t smearing = 0.005)
2913 {
2914   loadlibs();
2915
2916   TFile::Open("analysis_mc.root");
2917   hist = (TH1*) gFile->Get("dndeta_check_pt");
2918   
2919   dPhiHist = new TH1F("dPhiHist", ";#Delta phi", 400, -0.1, 0.1);
2920   dPhiHist2 = new TH1F("dPhiHist2", ";#Delta phi", 400, -0.1, 0.1);
2921   
2922   for (Int_t i=1; i<=hist->GetNbinsX(); i++)
2923   {
2924     Float_t pt = hist->GetBinCenter(i);
2925     Float_t deltaPhi = (0.076 - 0.039) / 2 / (pt / 0.15);
2926     
2927     if (smearing > 0)
2928     {
2929       gaus = new TF1("mygaus", "gaus(0)", -0.1, 0.1);
2930       gaus->SetParameters(1, -deltaPhi, smearing);
2931     
2932       dPhiHist->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
2933     
2934       dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
2935       gaus->SetParameters(1, deltaPhi, smearing);
2936       dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
2937     }
2938     else
2939 {
2940 dPhiHist->Fill(deltaPhi, hist->GetBinContent(i) / 2);
2941 dPhiHist2->Fill(deltaPhi, hist->GetBinContent(i) / 2);
2942 dPhiHist2->Fill(-deltaPhi, hist->GetBinContent(i) / 2);
2943 }
2944   }
2945   
2946   new TCanvas;
2947   dPhiHist->Draw();
2948   dPhiHist2->SetLineColor(2);
2949   dPhiHist2->Draw("SAME");
2950   gPad->SetLogy();
2951   
2952   TFile::Open("trackletsDePhi.root");
2953   //TFile::Open("tmp/correction_maponly-positive.root");
2954   //TFile::Open("tmp/correction_map.root");
2955   //tracklets = (TH1*) gFile->Get(Form("fDeltaPhi_%d", 0));
2956   tracklets = (TH1*) gFile->Get("DePhiPPTracklets");
2957   tracklets->Scale(1.0 / tracklets->GetMaximum() * dPhiHist->GetMaximum());
2958   tracklets->SetLineColor(4);
2959   tracklets->Draw("SAME");
2960 }
2961
2962 void VertexDistributions()
2963 {
2964   loadlibs();
2965   
2966   TFile::Open("correction_map.root");
2967   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
2968   if (!dNdEtaCorrection->LoadHistograms())
2969     return;
2970   
2971   all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all");
2972   trigger = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("trigger");
2973   vtx = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("vtx");
2974  
2975   nottriggered = (TH1*) all->Clone("nottriggered");
2976   nottriggered->Add(trigger, -1);
2977
2978   novertex = (TH1*) trigger->Clone("novertex");
2979   novertex->Add(vtx, -1);
2980   
2981   temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
2982   highmult = temphist->ProjectionX("highmult", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
2983   //all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
2984  
2985   for (Int_t i=1; i<=trigger->GetNbinsX(); i++)
2986   {
2987     all->SetBinContent(i, all->GetBinContent(i) / all->GetBinWidth(i));
2988     trigger->SetBinContent(i, trigger->GetBinContent(i) / trigger->GetBinWidth(i));
2989     vtx->SetBinContent(i, vtx->GetBinContent(i) / vtx->GetBinWidth(i));
2990     nottriggered->SetBinContent(i, nottriggered->GetBinContent(i) / nottriggered->GetBinWidth(i));
2991     novertex->SetBinContent(i, novertex->GetBinContent(i) / novertex->GetBinWidth(i));
2992     highmult->SetBinContent(i, highmult->GetBinContent(i) / highmult->GetBinWidth(i));
2993   }
2994
2995   new TCanvas;
2996   vtx->SetTitle("");
2997   vtx->SetStats(0);
2998   vtx->DrawCopy("HIST");
2999
3000   all->Scale(1.0 / all->Integral());
3001   nottriggered->Scale(1.0 / nottriggered->Integral());
3002   novertex->Scale(1.0 / novertex->Integral());
3003   highmult->Scale(1.0 / highmult->Integral());
3004
3005   new TCanvas;
3006   all->Draw("HIST");
3007   novertex->SetLineColor(2);
3008   novertex->Draw("HISTSAME");
3009   highmult->SetLineColor(3);
3010   highmult->Draw("HISTSAME");
3011
3012   legend = new TLegend(0.5, 0.5, 0.8, 0.8);
3013   legend->SetFillColor(0);
3014   legend->AddEntry(all, "all");
3015   legend->AddEntry(novertex, "no vertex");
3016   legend->AddEntry(highmult, "mult > 10");
3017   legend->Draw();
3018   
3019   new TCanvas;
3020   trigger->Scale(1.0 / trigger->Integral());
3021   vtx->Scale(1.0 / vtx->Integral());
3022   
3023   trigger->Divide(vtx);
3024   
3025   trigger->Draw();
3026   //vtx->SetLineColor(2);
3027   //vtx->Draw("SAME");
3028
3029   //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
3030   temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram();
3031   //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram();
3032   
3033   new TCanvas;
3034   legend = new TLegend(0.7, 0.7, 0.99, 0.99);
3035   legend->SetFillColor(0);
3036  
3037   Bool_t first = kTRUE; 
3038   for (Int_t i=0; i<20; i+=5)
3039   {
3040     highmult = temphist->ProjectionX("highmult", i+1, i+1+4);
3041     Printf("%f", highmult->Integral());
3042     if (highmult->Integral() <= 0)
3043       continue;
3044   
3045     for (Int_t j=1; j<=trigger->GetNbinsX(); j++)
3046       highmult->SetBinContent(j, highmult->GetBinContent(j) / highmult->GetBinWidth(j));
3047
3048     highmult->Scale(1.0 / highmult->Integral());
3049     highmult->SetLineColor((i/5)+1);
3050     highmult->GetYaxis()->SetRangeUser(0, 0.15);
3051     if (first)
3052     {
3053       highmult->DrawCopy();
3054       first = kFALSE;
3055     }
3056     else
3057       highmult->DrawCopy("SAME");
3058     legend->AddEntry(highmult->Clone(), Form("%d <= N <= %d", i, i+4));
3059   }
3060   legend->Draw();
3061  
3062 }
3063
3064 void PlotPt1DCorrection()
3065 {
3066   const char* files[] = { "field.root", "field_onlyprim.root", "nofield.root", "nofield_onlyprim.root" };
3067   const char* names[] = { "B: all", "B: primaries", "No B: all", "No B: primaries" };
3068   Int_t colors[] = { 1, 2, 3, 4 };
3069   
3070   loadlibs();
3071   
3072   dummy = new TH2F("dummy", ";p_{T};correction", 100, 0, 1.4, 100, 0.5, 3);
3073   dummy->SetStats(0);
3074   //dummy->GetYaxis()->SetTitleOffset(1.3);
3075   dummy->Draw();
3076   
3077   legend = new TLegend(0.48, 0.57, 0.88, 0.88);
3078   legend->SetFillColor(0);
3079   
3080   for (Int_t i=0; i<4; i++)
3081   {
3082     TFile::Open(files[i]);
3083     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
3084     if (!dNdEtaCorrection->LoadHistograms())
3085       return;
3086       
3087     hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get1DCorrectionHistogram("z", -9.9, 9.9, -0.79, 0.79);
3088     hist->SetLineColor(colors[i]);
3089     hist->SetLineWidth(2);
3090     hist->SetMarkerColor(colors[i]);
3091     hist->Draw("SAME");
3092     
3093     legend->AddEntry(hist, names[i], "L");
3094   }
3095   
3096   legend->Draw();
3097 }