0f9aede4ff28741ba7a404e6753a669e61e76b58
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / drawSystematicsNew.C
1
2 Int_t markers[] = {20,21,22,23,28,29};
3 Int_t colors[]  = {1,2,3,4,6,8,102};
4
5 void loadlibs()
6 {
7   gSystem->Load("libTree");
8   gSystem->Load("libVMC");
9
10   gSystem->Load("libSTEERBase");
11   gSystem->Load("libANALYSIS");
12   gSystem->Load("libPWG0base");
13 }
14
15 void InitPad()
16 {
17   if (!gPad)
18     return;
19
20   gPad->Range(0, 0, 1, 1);
21   gPad->SetLeftMargin(0.15);
22   //gPad->SetRightMargin(0.05);
23   //gPad->SetTopMargin(0.13);
24   //gPad->SetBottomMargin(0.1);
25
26   gPad->SetGridx();
27   gPad->SetGridy();
28 }
29
30 void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
31 {
32   //gROOT->ProcessLine(".L drawPlots.C");
33   gSystem->Load("libPWG0base");
34
35   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
36   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "dndeta_correction" };
37   const char* legendNames[] = { "#pi", "K", "p", "standard" };
38   Int_t folderCount = 3;
39
40   TH2F* h2DCorrections[4];
41   TH1F* h1DCorrections[4];
42   for (Int_t i=0; i<4; i++) {
43     TFile::Open(fileNames[i]);
44     AlidNdEtaCorrection* correctionTmp = new AlidNdEtaCorrection(folderNames[i],folderNames[i]);
45     correctionTmp->LoadHistograms();
46     
47     //    h2DCorrections[i] = correctionTmp->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get2DCorrectionHistogram("yz",-1,1);
48
49     h1DCorrections[i] = correctionTmp->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get1DCorrectionHistogram("z",-10,10,0.8,0.8);
50   }
51
52   TH2F* null = new TH2F("","",100,0.1,10,100,0.5,9.99);
53   null->SetXTitle("p_{T} (GeV/c)");
54   null->SetXTitle("Correction");
55
56   null->Draw();
57
58   //h1DCorrections[0]->SetMaximum(5);
59   //h1DCorrections[0]->Draw();
60   //  h2DCorrections[0]->Draw("colz");
61   for (Int_t i=1; i<4; i++) {
62     h1DCorrections[i]->Draw("same");
63   }
64   
65 }
66
67
68 void DrawEffectOfChangeInCrossSection(const char* fileName) {
69   
70   //get the data
71   TFile* fin = TFile::Open(fileName);
72
73   //const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore" };
74   const Char_t* changes[]  = { "pythia", "qgsm", "phojet" };
75   const Int_t nChanges = 3;
76   //const Char_t* changes[]  = {"pythia","ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25" };
77   Int_t colors[] = {1,1,kRed,kBlue,kGreen,kPink,1,kRed,kBlue};
78
79   TH1F* hRatios[3];
80   for(Int_t i=0; i<nChanges; i++) {
81     hRatios[i] = (TH1F*)fin->Get(Form("ratio_vertexReco_triggerBias_%s",changes[i]));
82     hRatios[i]->SetLineWidth(2);
83     hRatios[i]->SetLineColor(colors[i]);
84     hRatios[i]->SetMarkerStyle(22);
85     hRatios[i]->SetMarkerSize(0.8);
86
87     Float_t average = hRatios[i]->Integral(hRatios[i]->FindBin(-1), hRatios[i]->FindBin(1)) / (hRatios[i]->FindBin(1) - hRatios[i]->FindBin(-1) + 1);
88     Printf("%s: %.2f %%" , hRatios[i]->GetTitle(), (average - 1) * 100);
89   }
90
91   TPad* p = DrawCanvasAndPad("syst_changeInXsection",700,400);
92   p->SetRightMargin(0.2);
93   p->SetLeftMargin(0.13);
94
95   TH2F* null = new TH2F("","",100,-1.05,1.05,100,0.9,1.1);
96   null->GetXaxis()->SetTitle("#eta");
97   null->GetYaxis()->SetTitle("Ratio pythia/modified cross-sections");
98   null->Draw();
99
100   TLatex* text[3];
101
102   for(Int_t i=1; i<nChanges; i++) {
103     hRatios[i]->Draw("same");
104     
105     TString str(hRatios[i]->GetTitle());
106     str.Remove(0,str.First("DD"));
107     str.Remove(str.First(")"),1);
108     text[i] = new TLatex(1.08,hRatios[i]->GetBinContent(hRatios[i]->FindBin(0.))-0.002,str.Data());
109     text[i]->SetTextColor(colors[i]);
110     text[i]->Draw();
111   }
112 }
113
114
115 void DrawEffectOfChangeInComposition() {
116   
117   //get the data
118   TFile* fin = TFile::Open("systematics_composition.root");
119
120   Int_t colors[] = {1,2,103,102,4,6,1};
121
122   TH1F* hRatios[6];
123
124   Float_t etaRange = 0.899;
125
126   for (Int_t i=0; i<6; i++) {
127     hRatios[i] = (TH1F*)fin->Get(Form("ratio_%d",i));
128
129     hRatios[i]->SetLineWidth(2);
130     hRatios[i]->SetLineColor(colors[i]);
131     hRatios[i]->SetMarkerStyle(22);
132     hRatios[i]->SetMarkerSize(0.8);
133
134     hRatios[i]->GetXaxis()->SetRangeUser(-etaRange, etaRange);
135   }
136
137   TPad* p = DrawCanvasAndPad("syst_changeOfComposition",700,400);
138   p->SetRightMargin(0.2);
139   p->SetLeftMargin(0.13);
140
141   TH2F* null = new TH2F("","",100,-1.05,1.05,100,0.97,1.03);
142   null->GetXaxis()->SetTitle("#eta");
143   null->GetYaxis()->SetTitle("Ratio pythia/modified composition");
144   null->Draw();
145
146   TLatex* text[6];
147
148   for(Int_t i=0; i<6; i++) {
149     hRatios[i]->Draw("same");
150     
151     TString str(hRatios[i]->GetTitle());
152     str.Remove(0,16);
153     text[i] = new TLatex(1.08,hRatios[i]->GetBinContent(hRatios[i]->FindBin(0))-0.002,str.Data());
154     text[i]->SetTextColor(colors[i]);
155     text[i]->SetTextSize(0.053);
156     
157     text[i]->Draw();
158   }
159
160
161 }
162
163 TPad* DrawCanvasAndPad(const Char_t* name, Int_t sizeX=600, Int_t sizeY=500) {
164
165   gStyle->SetOptStat(0);
166   gStyle->SetOptTitle(0);
167   gStyle->SetOptFit(0);
168
169   gStyle->SetTextSize(0.04);
170   gStyle->SetTitleSize(0.05,"xyz");
171   //gStyle->SetTitleFont(133, "xyz");
172   //gStyle->SetLabelFont(133, "xyz");
173   //gStyle->SetLabelSize(17, "xyz");
174   gStyle->SetLabelOffset(0.01, "xyz");
175
176   gStyle->SetTitleOffset(1.1, "y");
177   gStyle->SetTitleOffset(1.1, "x");
178   gStyle->SetEndErrorSize(0.0);
179
180   //##############################################
181
182   //making canvas and pads
183   TCanvas *c = new TCanvas(name,name,sizeX,sizeY);
184
185   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
186
187   p1->SetBottomMargin(0.15);
188   p1->SetTopMargin(0.03);
189   p1->SetLeftMargin(0.15);
190   p1->SetRightMargin(0.03);
191   
192   p1->SetGridx();
193   p1->SetGridy();
194
195   p1->Draw();
196   p1->cd();
197
198   return p1;
199 }
200
201 void MisalignmentShowRawTrackPlots(const char* dirName = "fdNdEtaAnalysisESD")
202 {
203   loadlibs();
204
205   TFile* file = TFile::Open("fullA-simrec/MB2/analysis_esd_raw.root");
206   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
207   fdNdEtaAnalysis->LoadHistograms();
208
209   TFile* file2 = TFile::Open("fullA-sim/analysis_esd_raw.root");
210   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
211   fdNdEtaAnalysis2->LoadHistograms();
212
213   TH3* track1 = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track1");
214   TH3* track2 = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track2");
215
216   // normalize to number of events;
217   TH2* event1 = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
218   TH2* event2 = fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram();
219   Int_t event1Count = event1->Integral();
220   Int_t event2Count = event2->Integral();
221   track1->Scale(1.0 / event1Count);
222   track2->Scale(1.0 / event2Count);
223
224   const Float_t innerLimit = 0.49;
225   const Float_t outerLimit = 0.99;
226
227   track1->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
228   track2->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
229   AliPWG0Helper::CreateDividedProjections(track1, track2, "ze1");
230   TH1* fullRange = gROOT->FindObject("track1_ze1_div_track2_ze1");
231
232   track1->GetYaxis()->SetRangeUser(-innerLimit, innerLimit);
233   track2->GetYaxis()->SetRangeUser(-innerLimit, innerLimit);
234   AliPWG0Helper::CreateDividedProjections(track1, track2, "ze2");
235   TH1* central = gROOT->FindObject("track1_ze2_div_track2_ze2");
236   central->SetLineColor(1);
237   central->SetMarkerStyle(21);
238
239   for (Int_t x=1; x<track1->GetXaxis()->GetNbins(); ++x)
240     for (Int_t y=track1->GetYaxis()->FindBin(-innerLimit); y<track1->GetYaxis()->FindBin(innerLimit); ++y)
241       for (Int_t z=1; z<track1->GetZaxis()->GetNbins(); ++z)
242       {
243         track1->SetBinContent(x, y, z, 0);
244         track1->SetBinError(x, y, z, 0);
245         track2->SetBinContent(x, y, z, 0);
246         track2->SetBinError(x, y, z, 0);
247       }
248
249   track1->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
250   track2->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
251   AliPWG0Helper::CreateDividedProjections(track1, track2, "ze3");
252   TH1* peripheral = gROOT->FindObject("track1_ze3_div_track2_ze3");
253   peripheral->SetLineColor(2);
254   peripheral->SetMarkerStyle(22);
255   peripheral->SetMarkerColor(2);
256
257   TH2* tmp = new TH2F("tmp", ";p_{T} [GeV/c]      ;#frac{tracks full misalignment during rec.}{tracks ideal misalignment during rec.}", 1, 0.1, 10, 1, 0.9, 1.3);
258
259   tmp->SetStats(kFALSE);
260   //tmp->GetXaxis()->SetNoExponent();
261
262   Float_t ptStart = 0.1;
263
264   fullRange->GetXaxis()->SetRangeUser(ptStart, 9.9);
265   central->GetXaxis()->SetRangeUser(ptStart, 9.9);
266   peripheral->GetXaxis()->SetRangeUser(ptStart, 9.9);
267
268   TCanvas* canvas = new TCanvas("MisalignmentShowRawTrackPlots", "MisalignmentShowRawTrackPlots", 700, 400);
269   gPad->SetLogx();
270   gPad->SetGridx();
271   gPad->SetGridy();
272
273   TLegend* legend = new TLegend(0.2, 0.7, 0.4, 0.8);
274
275   legend->AddEntry(central, "|#eta| < 0.5");
276   legend->AddEntry(peripheral, "0.5 < |#eta| < 1.0");
277
278   legend->SetFillColor(0);
279
280   tmp->Draw();
281   //fullRange->Draw("SAME");
282   central->Draw("SAME");
283   peripheral->Draw("SAME");
284
285   legend->Draw();
286
287   canvas->SaveAs("syst_mis_ntracks.eps");
288 }
289
290
291 void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files, const char** dirs, const char** names, Int_t* histID)
292 {
293   loadlibs();
294   gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C");
295
296   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
297   canvas->Divide(2, 1);
298   canvas->cd(2)->SetGridx();
299   canvas->cd(2)->SetGridy();
300
301   TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
302   legend->SetFillColor(0);
303
304   TH1* base = 0;
305
306   for (Int_t i = 0; i < n; ++i)
307   {
308     TFile::Open(files[i]);
309
310     dNdEtaAnalysis* tmp = new dNdEtaAnalysis(dirs[i], dirs[i]);
311     tmp->LoadHistograms();
312
313     TH1* hist = tmp->GetdNdEtaPtCutOffCorrectedHistogram(histID[i]);
314
315     if (i == 0)
316       base = hist;
317
318     legend->AddEntry(hist, names[i]);
319
320     hist->SetMarkerColor(colors[i]);
321     hist->SetMarkerStyle(markers[i]);
322
323     canvas->cd(1);
324     hist->DrawCopy((i == 0) ? "" : "SAME");
325
326     if (i != 0)
327     {
328       PrintIntegratedDeviation(hist, base, names[i]);
329
330       canvas->cd(2);
331       hist->Divide(hist, base, 1, 1);
332       hist->GetYaxis()->SetRangeUser(0.9, 1.1);
333       hist->DrawCopy((i == 1) ? "" : "SAME");
334     }
335   }
336
337   legend->Draw();
338
339   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
340 }
341
342 void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files1, const char** files2, const char** dirs1, const char** dirs2, const char** names, Int_t* histID)
343 {
344   gSystem->Load("libPWG0base");
345
346   TCanvas* canvas = new TCanvas(canvasName, canvasName, 700, 400);
347   canvas->SetLeftMargin(0.12);
348
349   TLegend* legend = new TLegend(0.35, 0.7, 0.65, 0.85);
350   legend->SetFillColor(0);
351
352   for (Int_t i = 0; i < n; ++i)
353   {
354     TFile::Open(files1[i]);
355
356     dNdEtaAnalysis* tmp1 = new dNdEtaAnalysis(dirs1[i], dirs1[i]);
357     tmp1->LoadHistograms();
358
359     TFile::Open(files2[i]);
360
361     dNdEtaAnalysis* tmp2 = new dNdEtaAnalysis(dirs2[i], dirs2[i]);
362     tmp2->LoadHistograms();
363
364     TH1* hist1 = tmp1->GetdNdEtaPtCutOffCorrectedHistogram(histID[i]);
365     TH1* hist2 = tmp2->GetdNdEtaPtCutOffCorrectedHistogram(histID[i]);
366
367     TH1* division = hist1->Clone();
368
369     division->Divide(hist1, hist2, 1, 1, "B");
370
371     division->SetMarkerColor(colors[i]);
372     division->SetMarkerStyle(markers[i]);
373
374     legend->AddEntry(division, names[i]);
375
376     division->SetTitle("");
377     division->GetYaxis()->SetTitle("#frac{dN_{ch}/d#eta using MC vtx}{dN_{ch}/d#eta using ESD vtx}");
378     division->SetStats(kFALSE);
379     division->GetYaxis()->SetTitleOffset(1.3);
380     division->GetXaxis()->SetRangeUser(-0.99, 0.99);
381     division->GetYaxis()->SetRangeUser(0.981, 1.02);
382     division->DrawCopy((i == 0) ? "" : "SAME");
383   }
384
385   gPad->SetGridx();
386   gPad->SetGridy();
387
388   legend->Draw();
389
390   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
391 }
392
393 void MaterialBudgetChange()
394 {
395   const char* files[] =
396     { "Material-normal-mcvtx/analysis_esd.root",
397       "Material-increased-mcvtx/analysis_esd.root",
398       "Material-decreased-mcvtx/analysis_esd.root" };
399
400   const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
401   const char* names[] = { "no change", "+ 10 % material", "- 10 % material" };
402   Int_t hist[] = { 0, 0, 0 };
403
404   drawdNdEtaRatios("MaterialBudgetChange", 3, files, dirs, names, hist);
405 }
406
407 void MisalignmentChange()
408 {
409   const char* files[] =
410     { "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-sim-mcvtx/analysis_esd.root",
411       "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-simrec-mcvtx/analysis_esd.root" };
412
413   const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
414   const char* names[] = { "no change", "increased in sim/rec", "increased in sim" };
415   Int_t hist[] = { 0, 0, 0 };
416
417   drawdNdEtaRatios("MisalignmentChange", 2, files, dirs, names, hist);
418 }
419
420 void dNdEtaVertexRanges()
421 {
422   const char* files[] =
423     { "analysis_esd.root",
424       "analysis_esd.root",
425       "analysis_esd.root" };
426
427   const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
428   const char* names[] = { "|vtx-z| < 10 cm", "-10 cm < vtx-z < 0 cm", "0 cm < vtx-z < 10 cm" };
429   Int_t hist[] = { 0, 1, 2 };
430
431   drawdNdEtaRatios("dNdEtaVertexRanges", 3, files, dirs, names, hist);
432 }
433
434 void vertexShiftStudy(Int_t histID)
435 {
436   const char* files[] = { "maps/idealA/mc-vertex/analysis_esd.root", "results/idealC-idealA/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.05/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.1/analysis_esd.root", "maps/idealA/mc-vertex-shift-dep/analysis_esd.root" };
437   const char* dirs[] = { "dndeta", "dndeta", "dndeta", "dndeta", "dndeta" };
438   const char* names[] = { "mc vtx", "esd vtx", "+ 0.05 cm", "+ 0.1 cm", "old vtx shift" };
439   Int_t hist[] = { histID, histID, histID, histID, histID };
440
441   drawdNdEtaRatios("syst_vertex_shift1", 5, files, dirs, names, hist);
442
443   const char* files1[] = { "maps/idealA/mc-vertex/analysis_esd.root", "maps/idealA/mc-vertex/analysis_esd.root", "maps/idealA/mc-vertex/analysis_esd.root" };
444   const char* files2[] = { "results/idealC-idealA/analysis_esd.root", "results/idealC-idealA/analysis_esd.root", "results/idealC-idealA/analysis_esd.root" };
445   const char* dirs1[] = { "dndeta", "dndeta", "dndeta"};
446   const char* names[] = { "|vtx-z| < 10 cm", "-10 cm < vtx-z < 0 cm", "0 cm < vtx-z < 10 cm" };
447   Int_t hist[] = { 0, 1, 2 };
448
449   drawdNdEtaRatios("syst_vertex_shift2", 3, files1, files2, dirs, dirs, names, hist);
450 }
451
452 void vertexShift()
453 {
454   TFile::Open("vertex.root");
455
456   TH2* hist = gFile->Get("fVertexCorr");
457   TProfile* prof = hist->ProfileX();
458
459   prof->SetStats(kFALSE);
460   prof->SetTitle(";MC vtx-z [cm];mean (ESD vtx-z - MC vtx-z) [cm]");
461   prof->GetYaxis()->SetTitleOffset(1.2);
462
463   prof->SetLineWidth(2);
464
465   TCanvas* canvas = new TCanvas("syst_vertex_shift", "syst_vertex_shift", 700, 400);
466
467   gPad->SetGridx();
468   gPad->SetGridy();
469
470   prof->Draw();
471
472   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
473 }
474
475 void CompareRawTrackPlots(const char* fileName1, const char* fileName2, Float_t ptCut = 0.0, Int_t multCut = 1)
476 {
477   loadlibs();
478
479   const char* dirName = "fdNdEtaAnalysisESD";
480
481   TFile* file = TFile::Open(fileName1);
482   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
483   fdNdEtaAnalysis->LoadHistograms();
484
485   TFile* file2 = TFile::Open(fileName2);
486   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
487   fdNdEtaAnalysis2->LoadHistograms();
488
489   TH3* track1 = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track1");
490   TH3* track2 = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track2");
491
492   TH2* event1 = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
493   TH2* event2 = fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram();
494   Int_t event1Count = event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), multCut, event1->GetNbinsY() + 1);
495   Int_t event2Count = event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), multCut, event1->GetNbinsY() + 1);
496
497   Float_t nTrack1 = track1->Integral(track1->GetXaxis()->FindBin(-9.9), track1->GetXaxis()->FindBin(9.9), track1->GetYaxis()->FindBin(-0.79), track1->GetYaxis()->FindBin(0.79), track1->GetZaxis()->FindBin(ptCut), track1->GetZaxis()->GetNbins());
498   Float_t nTrack2 = track2->Integral(track2->GetXaxis()->FindBin(-9.9), track2->GetXaxis()->FindBin(9.9), track2->GetYaxis()->FindBin(-0.79), track2->GetYaxis()->FindBin(0.79), track2->GetZaxis()->FindBin(ptCut), track2->GetZaxis()->GetNbins());
499
500   Printf("%d tracks in %d events in first sample; %d tracks in %d events in second sample", (Int_t) nTrack1, event1Count, (Int_t) nTrack2, event2Count);
501
502   // normalize to number of events;
503   nTrack1 /= event1Count;
504   nTrack2 /= event2Count;
505
506   Printf("There are %.2f tracks/event in the first sample and %.2f tracks/event in the second sample. %.2f %% difference (with pt cut at %.2f GeV/c)", nTrack1, nTrack2, 100.0 * (nTrack1 - nTrack2) / nTrack1, ptCut);
507
508   gROOT->cd();
509
510   AliPWG0Helper::CreateDividedProjections(track1, track2);
511
512   new TCanvas; gROOT->FindObject("track1_yx_div_track2_yx")->Draw("COLZ");
513   new TCanvas; gROOT->FindObject("track1_zx_div_track2_zx")->Draw("COLZ");
514   new TCanvas; gROOT->FindObject("track1_zy_div_track2_zy")->Draw("COLZ");
515
516   for (Int_t i=0; i<3; i++)
517   {
518     char c = 'x' + (char) i;
519
520     /*proj1 = track1->Project3D(Form("%ce2", c));
521     proj2 = track2->Project3D(Form("%ce2", c));
522     AliPWG0Helper::NormalizeToBinWidth(proj1);
523     AliPWG0Helper::NormalizeToBinWidth(proj2);
524
525     new TCanvas;
526     proj1->DrawCopy();
527     proj2->SetLineColor(2);
528     proj2->SetMarkerColor(2);
529     proj2->DrawCopy("SAME");*/
530
531     AliPWG0Helper::CreateDividedProjections(track1, track2, Form("%ce2", c));
532     TH1* pt = gROOT->FindObject(Form("track1_%ce2_div_track2_%ce2", c, c));
533     new TCanvas; pt->DrawCopy();
534     gPad->SetGridx(); gPad->SetGridy();
535   }
536
537   event1_x = event1->ProjectionX("event1_x");
538   event1_y = event1->ProjectionY("event1_y");
539   event2_x = event2->ProjectionX("event2_x");
540   event2_y = event2->ProjectionY("event2_y");
541
542   new TCanvas; event1_x->DrawCopy(); event2_x->SetLineColor(2); event2_x->DrawCopy("SAME"); 
543
544   event1_x->Divide(event2_x);
545   event1_y->Divide(event2_y);
546
547   new TCanvas; event1_x->DrawCopy();
548   new TCanvas; event1_y->DrawCopy();
549
550   event1->Divide(event2);
551   new TCanvas;
552   event1->Draw("COLZ");
553   event1->SetMinimum(0.5);
554   event1->SetMaximum(2);
555
556
557 }
558
559 void MagnitudeOfCorrection(const char* fileName, const char* dirName = "dndeta", Float_t ptCut = 0.3)
560 {
561   loadlibs();
562
563   TFile* file = TFile::Open(fileName);
564   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
565   fdNdEtaAnalysis->LoadHistograms();
566   fdNdEtaAnalysis->GetData()->PrintInfo(ptCut);
567 }
568
569 Double_t ConvSigma1To2D(Double_t sigma)
570 {
571   return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
572 }
573
574 Double_t ConvDistance1DTo2D(Double_t distance)
575 {
576   return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
577 }
578
579 Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
580 {
581   Double_t count = 0;
582
583   //nSigma = ConvSigma1To2D(nSigma);
584
585   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
586     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
587     {
588       Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
589       Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
590
591       Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
592
593       d = ConvDistance1DTo2D(d);
594
595       if (d < nSigma)
596         count += tracks->GetBinContent(x, y);
597     }
598
599   return count;
600 }
601
602 TH2F* Sigma2VertexGaussianTracksHist()
603 {
604   TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
605
606   TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
607   gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
608
609   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
610     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
611       tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
612
613   //normalize
614   tracks->Scale(1.0 / tracks->Integral());
615
616   return tracks;
617 }
618
619 TH1F* Sigma2VertexGaussian()
620 {
621   TH2F* tracks = Sigma2VertexGaussianTracksHist();
622
623   TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
624   canvas->Divide(2, 2);
625
626   canvas->cd(1);
627   tracks->Draw("COLZ");
628
629   TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 50, 0.05, 5.05);
630   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
631     ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
632   ratio->SetMarkerStyle(21);
633
634   canvas->cd(2);
635   ratio->DrawCopy("P");
636
637   TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 4 sigma / % included n sigma", 50, 0.05, 5.05);
638   Double_t sigma3 = Sigma2VertexCount(tracks, 4);
639   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
640     ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma)));
641   ratio2->SetMarkerStyle(21);
642
643   canvas->cd(3);
644   ratio2->DrawCopy("P");
645
646   canvas->SaveAs("Sigma2Vertex.eps");
647
648   return ratio2;
649 }
650
651 TH1F** Sigma2VertexSimulation(const char* fileName = "systematics.root")
652 {
653   TFile* file = TFile::Open(fileName);
654
655   TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertexTracks"));
656   TH1F* sigmavertexPrim = dynamic_cast<TH1F*> (file->Get("fSigmaVertexPrim"));
657   if (!sigmavertex || !sigmavertexPrim)
658   {
659     printf("Could not read histogram(s)\n");
660     return;
661   }
662
663   // calculate ratio
664   TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included in 4 #sigma / % included in N#sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
665
666   // calculate contamination
667   TH1F* contamination = ratio->Clone("sigmavertexsimulation_contamination");
668   contamination->SetTitle("sigmavertexsimulation_contamination;N#sigma;1 + N_{secondaries} / N_{all}");
669
670   for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
671   {
672     ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(4)) / sigmavertex->GetBinContent(i));
673     contamination->SetBinContent(i, 1 + (sigmavertex->GetBinContent(i) - sigmavertexPrim->GetBinContent(i)) / sigmavertex->GetBinContent(i));
674   }
675
676   // print stats
677   for (Float_t sigma = 2.0; sigma < 5.25; sigma += 0.5)
678   {
679     Float_t error1 = 1 - ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma - 0.5));
680     Float_t error2 = -1 + ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma + 0.5));
681     Float_t cont = -1 + contamination->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma));
682     Printf("%.2f sigma --> syst. error = - %.2f %% + %.2f %%, cont. = %.2f %%", sigma, error1 * 100, error2 * 100, cont * 100);
683   }
684
685   TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
686   canvas->Divide(2, 1);
687
688   canvas->cd(1);
689   sigmavertex->SetMarkerStyle(21);
690   sigmavertex->Draw("P");
691
692   canvas->cd(2);
693   ratio->SetMarkerStyle(21);
694   ratio->DrawCopy("P");
695
696   contamination->DrawCopy("SAME");
697
698   TH1F** returnContainer = new TH1F*[2];
699   returnContainer[0] = ratio;
700   returnContainer[1] = contamination;
701
702   return returnContainer;
703 }
704
705 void Sigma2VertexCompare(const char* fileName = "systematics.root")
706 {
707   TH1F* ratio1 = Sigma2VertexGaussian();
708
709   TH1F** hists = Sigma2VertexSimulation(fileName);
710   TH1F* ratio2 = hists[0];
711   TH1F* contamination = hists[1];
712
713   ratio1->SetStats(kFALSE);
714   ratio2->SetStats(kFALSE);
715
716   ratio1->SetMarkerStyle(0);
717   ratio2->SetMarkerStyle(0);
718
719   ratio1->SetLineWidth(2);
720   ratio2->SetLineWidth(2);
721
722   TLegend* legend = new TLegend(0.6, 0.8, 0.95, 0.95);
723   legend->SetFillColor(0);
724   legend->AddEntry(ratio1, "Gaussian");
725   legend->AddEntry(ratio2, "Simulation");
726   legend->AddEntry(contamination, "1 + Contamination");
727
728   ratio2->SetTitle("");
729   ratio2->GetYaxis()->SetTitleOffset(1.5);
730   ratio2->GetXaxis()->SetRangeUser(2, 5);
731
732   TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
733   InitPad();
734
735   ratio2->SetMarkerStyle(21);
736   ratio1->SetMarkerStyle(22);
737
738   ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
739   ratio2->SetLineColor(kRed);
740   ratio2->SetMarkerColor(kRed);
741   ratio2->Draw("PL");
742   ratio1->Draw("SAMEPL");
743
744   contamination->Draw("SAME");
745
746   legend->Draw();
747
748   canvas->SaveAs("Sigma2VertexCompare.eps");
749 }