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